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Abstract 

As  theory  lags  experiment  for  dielectric  barrier  discharge  flow  control,  two 
different  computational  methods  are  implemented  to  give  further  insight  into 
characteristics  of  the  dielectric  barrier  discharge  (DBD).  A  one-dimensional  fluid  model 
of  a  surface-type  dielectric  barrier  discharge  is  created  using  He  as  the  background  gas. 
This  simple  model,  which  only  considers  ionizing  collisions  and  recombination  in  the 
electropositive  gas,  creates  an  important  framework  for  future  studies  into  the  origin  of 
experimentally  observed  flow-control  effects  of  the  DBD.  The  two  methods  employed  in 
this  study  include  the  semi-implicit  sequential  algorithm  and  the  fully  implicit 
simultaneous  algorithm.  The  first  involves  consecutive  solutions  to  Poisson’s,  the 
electron  continuity,  ion  continuity  and  electron  energy  equations.  This  method  combines 
a  successive  over-relaxation  algorithm  as  a  Poisson  solver  with  the  Thomas  algorithm 
tridiagonal  routine  to  solve  each  of  the  continuity  equations.  The  second  algorithm 
solves  an  Ax=b  system  of  linearized  equations  simultaneously  and  implicitly.  The 
coefficient  matrix  for  the  simultaneous  method  is  constructed  using  a  Crank-Nicholson 
scheme  for  additional  stability  combined  with  the  Newton-Raphson  approach  to  address 
the  non-linearity  and  to  solve  the  system  of  equations.  Various  boundary  conditions,  flux 
representations  and  voltage  schemes  are  modeled.  Test  cases  include  modeling  a 
transient  sheath,  ambipolar  decay  and  a  radio-frequency  discharge.  Results  are  compared 
to  validated  computational  solutions  and/or  analytic  results  when  obtainable.  Finally,  the 
semi-implicit  method  is  used  to  model  a  DBD  streamer. 
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COMPUTATIONAL  MODEL  OF  ONE-DIMENSIONAL 


DIELECTRIC  BARRIER  DISCHARGES 


For  this  project,  the  interest  in  the  dielectric  barrier  discharge  (DBD)  stems  from 
its  application  to  the  plasma  actuator.  The  plasma  actuator  is  a  flow  control  device  with 
no  moving  parts,  does  not  change  airfoil  shape  and  puts  no  parts  in  the  flow.  Many 
experimental  endeavors  are  already  underway  studying  the  capabilities  of  the  plasma 
actuator  (1;2;3).  Observation  without  physical  explanation  leaves  much  to  be  desired. 
The  computational  side  of  this  technology  has  yet  to  fully  describe  the  plasma  actuator 
phenomena  theoretically.  This  project  tests  and  validates  a  numerical  model  that  is 
ultimately  used  to  simulate  a  one-dimensional  DBD  in  order  to  study  the  effect  of  the 
dielectric  barrier  on  the  discharge. 

I.  DBD  History  and  Applications 

The  DBD  has  an  extensive  background  in  both  industrial  and  scientific 
applications.  While  the  basic  theory  of  the  DBD  can  be  found  in  most  plasma  dynamics 
texts  (4;  5),  flow  control  DBD  devices  extend  operations  beyond  the  domain  of  current 
theory.  The  scientific  community  is  currently  attempting  to  remedy  this  shortcoming. 
This  chapter  gives  a  brief  history  of  the  DBD,  introduces  the  basic  theory  and  reviews 
recent  experimental  and  computational  endeavors. 
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Brief  History 

Much  of  the  research  on  “silent  discharges”  -  another  name  for  the  DBD  -  finds 
its  roots  in  the  evolution  of  industrial  processes.  The  original  experimental  set-up  driving 
the  DBD  can  be  traced  to  1857  when  Siemens  proposed  a  special  electrical  discharge  to 
produce  ozone  (6:309).  In  1955,  Tanaka  discovered  that  the  DBD  could  be  used  for 
excimer  fonnation  in  rare  gases  (6:309).  By  1960,  Bitzer  and  Slottow  applied  this 
discovery  to  the  invention  of  the  plasma  display  panel  (PDP)  (7:R54).  First  developed  as 
a  monochrome  display  for  educational  purposes,  the  PDP  has  nearly  superseded  the 
cathode  ray  tube  for  retail  color  television  technology.  While  current  DBD  investigations 
include  a  wide  variety  of  applications,  one  of  significant  interest  to  the  Air  Force  is  DBD 
flow  control. 

AF  Applications 

With  contemporary  Air  Force  weapons  systems  progressing  toward  smaller, 
faster,  smarter  designs,  any  scientific  endeavor  that  could  possibly  further  that  evolution 
is  on  the  forefront  of  the  military  scope.  Compared  with  the  now-archaic  BLU-109 
bomb,  the  smart  bomb  offers  identical  penetration  capabilities  in  a  package  one-third  the 
diameter,  half  the  length  and  only  an  eighth  of  the  total  weight  (8:7).  As  things  get 
smaller,  they  get  faster  and  smarter  simultaneously.  Current  operations  require  GPS 
equipped  small-diameter  bombs  with  top-of-the-line  control  systems  ranging  from  lattice 
folding  fins  to  laser  guidance  (8:7).  These  recent  developments  have  given  US  military 
forces  the  ability  to  do  things  never  dreamed  of  in  the  past. 
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As  more  and  more  warfighters  express  an  interest  in  continuing  to  enhance  the 
agility  of  their  weapons,  technology  simultaneously  approaches  its  miniaturization  limit 
for  control  systems.  A  new  flow  control  device  would  greatly  enhance  the  ability  to 
continue  decreasing  the  size  and  increasing  the  control  capabilities  of  future  combat 
technologies.  The  DBD  plasma  actuator  could  possibly  be  a  part  of  this  future.  While 
actuator  flow  control  is  a  very  advanced  subject,  the  physical  laws  driving  the  device 
originate  in  the  basic  physics  of  plasma  dynamics. 

Plasma  Dynamics 

In  general,  a  gas  discharge  is  generated  by  power  input  to  a  background  gas 
maintained  in  a  capacitive  set-up  between  two  electrodes  (Fig.  1).  The  voltage  source  can 
be  either  AC  or  DC.  For  each  source,  the  potential  difference  (or  alternating  potential 
difference)  between  the  electrodes  establishes  an  electric  field.  The  field  accelerates 
electrons  -  either  free  in  the  gas  or  pulled  from  the  cathode  -  through  the  background 
gas.  Some  of  these  translational  electrons  collide  with  the  neutral  particles  in  the 
background  gas  generating  ion/electron  pairs. 

Background  Gas 

@  0~V 

Figure  1 .  Basic  Plasma  Discharge  Configuration 

At  a  certain  value  of  applied  potential,  the  electrons  gain  enough  energy  to  create 
an  electron  avalanche.  One  seed  electron  accelerates  through  the  field  to  a  high  enough 
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energy  to  knock  an  additional  electron  off  neutral  particles  or  other  positive/negative 
ions.  These  additional  electrons  accelerate  and  cause  even  more  ionization.  This  process 
is  known  as  breakdown  and  the  potential  difference  required  to  initiate  it  is  the 
breakdown  voltage. 

A  self-sustaining  discharge  occurs  when  the  applied  voltage  is  large  enough  to 
initiate  breakdown.  In  any  type  of  self-sustaining  discharge,  the  cathode  sheath  plays  an 
integral  role  because  most  of  the  ionization  sustaining  the  discharge  occurs  in  this  region. 
The  size  of  the  sheath  is  very  small  and  depends  on  the  type  and  pressure  of  the 
background  gas  (9: 134).  The  cathode  sheath  of  a  DC  discharge  has  net  positive  charge 
that  essentially  shields  the  rest  of  the  discharge  from  the  cathode  potential  (10:9).  The 
effective  potential  difference  for  the  non-sheath  plasma  is  greatly  diminished  because  of 
the  significant  potential  drop  that  occurs  in  the  sheath. 


Figure  2.  Current/Voltage  Discharge  Classification  (11 :84) 
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Once  breakdown  has  occurred  and  the  sheath  has  been  established,  gas  discharges 
are  classified  according  to  characteristic  voltage,  current  and  pressure  ranges  as  shown  in 
Fig.  2.  Since  the  DBD  is  a  special  case  of  the  glow  discharge,  this  region  of  the  chart  of 
discharge  classification  is  of  most  interest  to  this  project.  The  gas  in  a  glow  discharge  is 
only  weakly  ionized  (10:9)  but  the  ionization  has  not  dropped  low  enough  to  lose  its 
plasma  characteristics.  The  common  qualities  of  low  pressures  ( — 1-10  Torr),  low 
currents  (~10"  -  10'  A)  and  high  voltages  (~1  (T-  1 0  V)  generally  distinguish  this  type 
of  discharge  (9:2-4). 

Dielectric  Barrier  Discharge 

The  DBD  finds  its  home  near  the  high-end  of  the  glow  discharge  pressure 
spectrum.  As  the  name  suggests,  a  dielectric  barrier  is  inserted  into  the  discharge  system 
by  coating  one  or  both  of  the  electrodes.  This  barrier  blocks  all  (or  most)  of  the  current 
flow  to  the  buried  electrode  and  alters  the  electric  potential  across  the  discharge. 


>AM - 1| — \ 

R plasma  Cdielectric 


V applied 

Figure  3.  Equivalent  Circuit  Model  of  DBD 

Effectively,  the  dielectric  barrier  creates  a  capacitance  in  the  discharge  circuit.  As 
an  example  in  Fig.  3,  if  a  large  negative  potential  is  applied  to  the  left  electrode  then 
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electrons  will  accelerate  toward  the  right  electrode.  If  this  electrode  is  covered  with  a 
dielectric  that  blocks  current  and  accumulates  this  charge,  this  scenario  is  analogous  to 
charging  a  capacitor  in  a  DC  circuit.  The  electron  current,  had  the  barrier  not  been 
inserted,  turns  into  displacement  current.  Assuming  that  the  plasma  has  some 
characteristic  resistance,  the  discharge  becomes  a  capacitive  circuit  (Fig.  3). 

The  charge  build-up  on  the  dielectric  surface  dramatically  changes  the  discharge 
behavior.  When  the  applied  voltage  is  DC,  the  charge  accumulation  eventually  reduces 
the  effective  electric  field  eventually  driving  the  field  to  zero  and  extinguishing  the 
discharge  all-together.  In  an  AC  discharge,  the  charge  accumulation  on  the  dielectric 
surface  may  give  a  large  number  of  seed  electrons  for  the  backward  stroke  (as  opposed  to 
the  assumed  few  seed  electrons  on  the  forward  stroke).  Some  numerical  models  predict 
that  this  asymmetry  is  the  source  of  the  experimentally  observed  flow  control  (12:9; 
13:11). 

Recent  Experimental  Work 

Experimental  research  has  proven  that  the  DBD  can  be  used  as  a  flow  control 
device.  In  wind  tunnels,  an  actuator-equipped  airfoil  has  been  shown  to  attain  higher 
angles  of  attack  before  separation  and  stall  occur  (1:3).  This  coincides  with  the  increased 
coefficient  of  lift  associated  with  having  the  actuator  turned  on  (Fig. 4).  In  this 
experiment,  a  higher  coefficient  of  drag  was  simultaneously  observed  when  only  one 
actuator  was  used.  However,  when  four  actuators  were  placed  in  series  along  the  camber, 
the  increased  coefficient  of  lift  was  maintained  while  the  increased  drag  was  essentially 
eliminated. 
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Experiment  has  also  documented  that  the  actuator  can  reattach  separated  flow 
(2:8).  Fig.  5  shows  the  flow  visualization  created  by  introduction  of  smoke  streaklines 
into  the  flow.  The  top  picture  shows  the  separated,  turbulent  flow  under  the  airfoil  and 
the  bottom  picture  shows  reattached,  laminar  flow.  The  only  configuration  change 
between  the  two  photos  is  the  actuator  power. 


NAGA  OOC9 


Figure  4.  Actuator-On  Fift  Increase  (1:3) 


Figure  5.  Actuator  Reattachment  of  Flow  (2:8) 
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(b) 


Figure  6.  Induced  Airflow  in  Initially  Still  Air  (3:2125) 


Upper 

Electrode 


Dielectric 

I  -I  Lower  "" 

Electrode 

Figure  7.  Asymmetric  Electrode 

Finally,  the  actuator  has  been  shown  to  induce  flow  in  initially  still  air.  Fig.  6 
shows  titanium  dioxide  smoke  emitted  from  a  vertical  tube  2.5mm  above  a  series  of 
asymmetric  electrodes  (3:2125).  Each  asymmetric  electrode  consists  of  one  exposed 
electrode  sitting  on  top  of  the  dielectric  and  one  electrode  sitting  immediately  adjacent  on 
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the  bottom  surface  of  the  dielectric  as  shown  in  Fig.  7.  Each  asymmetric  electrode 
induces  neutral  gas  flow.  With  the  electrodes  phased  so  that  they  can  accelerate  the  air 
flow,  the  smoke  moves  to  the  right  or  left  according  to  the  phasing  of  the  electrodes.  The 
only  configuration  change  between  the  two  pictures  in  Fig.  6  is  the  phasing. 

Experimental  results  definitively  show  that  the  DBD  plasma  actuator  influences 
or  creates  airflow.  In  order  to  maximize  actuator  capabilities,  however,  scientists  must  be 
able  to  theoretically  explain  the  results.  Numerous  computational  groups  are  working  to 
discover  the  source  and  nature  of  airflow  alteration  when  the  DBD  actuator  is  engaged. 

Numerical  Efforts 

Creating  a  numerical  model  of  DBD  flow  control  is  a  difficult  endeavor.  Keeping 
track  of  huge  particle  densities  (~1015-  102°  particles/cm3),  differing  species,  a  variety  of 
possible  particle  relations  (momentum  transfer,  excitations,  ionization,  electrical 
interactions)  and  varying  electric  field  configurations  proves  to  be  a  significant  challenge. 
In  the  past  twenty  years  with  the  advent  of  continually  improving  computers  and 
processors,  many  numerical  models  have  made  their  way  into  academic  periodicals  each 
building  on  and/or  fine-tuning  some  previous  set  of  calculations. 

While  there  are  far  too  many  models  to  cover  each  in  depth,  this  project  primarily 
utilizes  the  work  of  Boeuf  and  Pitchford,  Hilbun  and  Font  (14;  15;  16;  12)  as  guides  for 
the  numerical  model  and  validation  tests.  Two  Boltzmann  solvers,  BOLSIG  and  SIGLO- 
RF  (17;  18),  are  used  to  generate  the  rates/transport  coefficients  as  well  as  results  for 
qualitative  comparison.  BOLSIG  determines  the  solution  to  the  Boltzmann  equation  for 
electrons  in  weakly  ionized  gases.  These  solutions  assume  steady-state,  uniform  fields 
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and  can  be  calculated  for  fifteen  different  gases  and  a  wide  range  of  reduced  field  values 
(19:1).  SIGLO-RF  is  a  one-dimensional  fluid  model  of  a  radio-frequency,  low  pressure 
discharge.  This  code  gives  the  ion  density,  electron  density  and  electric  field  as  functions 
of  time  and  one-dimensional  space  (20: 1).  Once  again,  fifteen  different  gases  can  be 
modeled. 

So  far,  computational  projects  have  only  calculated  a  relatively  small  asymmetry 
associated  with  the  DBD.  Font  reports  a  time  averaged  computational  force  of  1.5  x  1 0 
N  (12:9)  for  an  actuator  250mm  wide.  This  force  is  roughly  equivalent  to  the  weight  of 
l/40th  of  a  teaspoon  of  salt.  The  importance  of  the  DBD  model  greatly  increases 
knowing  that  the  search  could  be  likened  to  measuring  the  effect  of  several  grains  of  salt 
on  the  wind.  Using  basic  fluid  dynamic  principles  to  model  an  actuator  in  one- 
dimension,  exploring  both  AC  and  DC  field  sources  and  the  addition  of  a  dielectric 
substrate,  this  project  will  help  further  establish  the  dielectric  barrier  effect  on  discharge 
current.  This  paper  covers  the  basic  theory  and  equations,  computational  development, 
numerical  validation  and  results  for  the  DBD  streamer  simulation. 
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II.  Fluid  Approach  to  DBD  Simulation 


While  there  are  several  methods  that  could  be  employed  to  model  the  DBD,  this 
project  utilizes  the  fluid  approach.  The  fluid  approach  assumes  that  microscopic  detail  is 
unnecessary  for  determining  macroscopic  behavior.  The  background  theory  and  model 
development  are  covered  in  this  chapter.  Unless  otherwise  specified,  all  variables  and 
equations  are  expressed  in  tenns  of  SI  units. 


The  Boltzmann  Equation 

Like  many  problems  relating  to  plasmas,  this  story  begins  with  the  Boltzmann 
equation 


dA 

dt 


+  v-Vfs+a-V,fs 


dt 


J  collission 


(1) 


The  solution  to  this  equation,  fs  (x,  y,  z,  vx ,  v  ,vz,t),  gives  the  time-dependent  single 

particle  distribution  function  in  phase  space  (21:11).  A  particle  in  Eq.  ( 1)  is  of  type  “s,” 
has  velocity  v  and  acceleration  a  .  The  subscript  “s”  could  designate  electrons,  positive 
ions,  negative  ions  or  neutrals.  The  first  gradient,  V/s ,  represents  the  change  to  the 
distribution  function  with  respect  to  configuration  coordinates,  (x,y,z),  whereas  the 
second  gradient,  V  vfs ,  represents  the  change  to  the  distribution  function  with  respect  to 
the  velocity  coordinates,  (vx,vv,v_)  .  The  right-hand  side  accounts  for  the  changes  to  the 
distribution  function  brought  about  by  collisions  (2 1:12). 
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The  solution  to  the  Boltzmann  equation  represents  the  number  of  particles  in  a 
volume  element  dx-dy  ■  dz  at  position  x  =  (x,  v,  z)  and  the  velocity  volume  element 
dvx  ■  dv  ■  dv,  with  velocity  v  =  (vx,v  ,vz)  at  time  t  (22:30).  This  solution  gives  a 

detailed,  microscopic  view  of  any  given  plasma  system  at  any  given  time.  Many  times, 
due  to  the  complexity  of  associated  systems,  it  is  impossible  to  find  an  analytic  solution 
and/or  unwieldy  to  calculate  for  every  time  step  in  the  process.  When  the  system  under 
investigation  does  not  require  microscopic  detail,  the  fluid  equations  provide  a  more 
practical  alternative  for  solution. 


The  Fluid  Equations 

The  fluid  equations  are  derived  by  taking  velocity  moments  of  the  Boltzmann 
equation.  The  nth  velocity  moment  can  be  represented  as  (2 1:19): 


^-  +  v-V/,+o-V,/,= 


(§2 

V  St  J 

collission  _ 

dvxdvydv. 


(2) 


These  moments  generate  equations  in  terms  of  macroscopic  variables  that  can  be 
measured  or,  in  the  case  of  this  study,  solved  for  computationally.  The  details  of  their 
derivation  will  not  be  covered  here  but  can  be  found  in  elementary  plasma  references 
(21:19-26;  22:39-42). 

The  model  assumptions  become  important  to  the  simplification  of  each  of  the 
moments.  The  model  assumptions  for  this  project  along  with  their  associated 
simplifications  are  discussed  in  Appendix  A.  The  simplified  zeroth  moment  yields  the 
particle  continuity  equation: 
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(3) 


where  ns  is  the  density  of  the  particle  of  type  “s,”  vs  is  average  particle  velocity  and  Ss 

is  the  particle  source/loss  function  (22:39).  The  first  moment  gives  the  particle 
momentum  equation,  which  can  be  simplified  to  yield  the  drift-diffusion  flux 
approximation: 

f)  =  nsvs  =  +njisE  -  DsVns  (4) 

where  Ts  represents  the  flux  of  “s”  type  particles,  E  is  the  electric  field,  +  jus  is  the 
mobility  coefficient  and  Ds  is  the  free  diffusion  coefficient  of  particle  of  type  “s” 

(14: 1378).  The  selected  sign  of  the  mobility  tenn  matches  that  of  the  charge  of  the 
particle.  Finally,  the  second  moment  is  the  electron  energy  equation  which  in  simplified 
form  greatly  resembles  Eq.  (3): 

nJfA  +  ly.f£=S£  (5) 

8t  3 

where  ue  represents  the  average  electron  energy,  Se  represents  the  electron  energy 
source/loss  function  (10:4).  The  electron  energy  density  flux  is  given  by 

f e  =  ~{neue  )neE  -  DeV{neue ) .  (6) 

The  continuity  and  momentum  equations  are  considered  sufficient  to  model  ion 
transport  because  the  local  field  approximation  is  applied  (10:14).  This  approximation 
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assumes  a  direct  relationship  between  the  particle  energy  distribution  and  the  electric 
field.  For  the  electrons,  however,  this  approximation  is  questionable.  Therefore,  the 
electron  energy  can  not  be  directly  related  to  the  local  field  and  Eq.  (5)  must  be  solved  in 
addition  to  Eqs.  (3)  and  (4)  in  order  to  more  accurately  characterize  electron  behavior. 

The  computational  variables  for  this  project  are  ne,  Te,  n  ,  f  ,  ( neue ),  Te  and 

E .  While  Eqs.  (3)-(6)  give  us  the  ability  to  solve  for  the  first  six  of  these  variables,  the 
transport  parameters,  source/loss  terms  and  vital  sixth  variable  remain  undetermined. 
Evaluation  of  the  electric  field  requires  a  self-consistent  solution  of  Poisson’s  equation 
and  will  be  covered  in  the  next  section.  The  transport  coefficients  ( jus  and  Ds )  and  the 

production/loss  rates,  essential  to  an  accurate  form  of  the  fluid  equations,  will  be  covered 
immediately.  Note  that  for  this  study,  the  values  for  all  transport  coefficients  and 
source/loss  rates  are  evaluated  by  a  numerical  solution  of  the  collisional  Boltzmann 
equation  (17;  18). 

Transport  Coefficients 

Transport  coefficients  specify  the  macroscopic  properties  of  the  plasma  and  are 
determined  from  a  selected  weighting  of  the  distribution  function.  These  coefficients  are 
typically  expressed  as  functions  of  average  energy  or  the  reduced  electric  field  expressed 
as  E/N  or  El  p  (9:17).  The  reduced  field  in  terms  of  the  number  density  ( E  /  N  )  has 
units  of  Townsend  (Td)  where 


l-Td  =  \x\0  11 V  •  cm2 .  (7) 

Ignoring  local  pressure  variations,  E  /  77  is  related  to  E  /  p  by  the  ideal  gas  law 
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p  =  NkBT 


(8) 


where  p  is  the  pressure  of  the  gas,  N  is  the  neutral  number  density,  kB  is  Boltzmann’s 
constant  and  T  is  the  temperature  of  the  gas.  Assuming  the  temperature  is  3 00°  A"  then 

3  •  Td  « 1  •  F  /  cm  /  torr  .  (9) 

The  transport  coefficients  each  apply  to  one  term  of  the  drift-diffusion  flux 
approximation: 

rs  =  +  nsMsE  -  •  (10) 

1  2 

In  the  first  tenn  of  Eq.  (10),  the  mobility  coefficient  jus  relates  the  proportionality 
between  the  drift  velocity  of  a  charged  particle  and  the  field: 

Vd,e=-Me[Ue]E  O1) 

and 

^d,P  =  Mp[E/ p]E  (12) 

where  vd  e  is  the  drift  velocity  of  the  electron  (9: 11),  vd  p  is  the  drift  velocity  of  the  ion, 
pe  [ue  ]  is  the  electron  mobility  based  on  the  local  average  electron  energy,  E  /  p  is  the 
local  reduced  field  expressed  in  terms  of  pressure  and  pp  [ E  /  p\  is  the  ion  mobility  based 
on  the  local  reduced  field. 
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Figure  8.  Electron  Mobility  Coefficient  in  Helium  at  1  Torr  (17) 


Figure  9.  Ion  Mobility  Coefficient  in  Helium  at  1  Torr  (18) 


The  difference  between  the  ion  and  electron  mobility  functional  dependence  is 
once  again  based  on  the  local  field  approximation.  Because  the  ion  energy  distribution 
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can  be  directly  related  to  the  local  field,  the  ion  mobility  is  considered  to  be  a  function  of 
the  reduced  field.  The  electron  mobility  is  a  function  of  the  average  electron  energy 
because  the  local  field  approximation  does  not  apply. 

Figs.  8  and  9  give  examples  of  the  ion  and  electron  mobility  in  helium  gas.  The 
plot  domains  are  limited  to  the  energy  and  reduced  field  values  that  are  most  commonly 
encountered  in  a  helium  discharge.  Also  note  from  Figs.  8  and  9  that  the  electron 
mobility  is  two  to  three  orders  of  magnitude  greater  than  the  ion  mobility. 

The  magnitude  of  the  electron  mobility  contributes  significantly  to  the  stiffness  of 
the  electron  equations  because  the  response  of  the  electron  velocity  to  the  field  can 
change  a  great  deal  over  a  small  spatial  range.  The  reduced  magnitude  of  the  ion 
mobility  means  that  the  ions  are  not  nearly  as  responsive  to  the  field.  Because  of  this, 
many  numerical  models  use  either  a  constant  ion  mobility  or  an  equation  that 
characterizes  the  ion  mobility  for  certain  large  ranges  of  the  reduced  field  (15:561 1; 
23:2789). 

The  transport  coefficient  associated  with  part  2  of  Eq.  (10)  is  the  diffusion 
coefficient  D .  The  diffusion  coefficient  is  determined  by 


D  = 


(13) 


where  v  is  the  velocity  of  the  particle  and  vm  is  the  effective  collision  frequency  for 

momentum  transfer  (5:9;  20).  This  parameter  quantifies  the  flux  associated  with  the 
spatial  non-uniformity  of  the  particle  or  energy  density  in  the  gas.  The  electron  diffusion 
coefficient  in  helium  is  shown  in  Fig.  10. 
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Figure  10.  Electron  Diffusion  Coefficient  in  Helium  at  1  Torr  (17) 


The  diffusion  coefficient  for  electrons  is  assumed  to  be  a  function  of  the  average 
energy  determined  by  the  energy  continuity  equation.  The  ion  diffusion  coefficient  is  is  a 
function  of  E  /  p  and  directly  related  to  the  ion  mobility  coefficient.  Since  the  ions  are 
assumed  to  have  a  Maxwellian  distribution,  the  Einstein  relation  can  be  employed: 


D  _  kBT 
/u  e 


(14) 


where  e  is  the  unsigned  electron  charge,  kB  is  Boltzmann’s  constant,  T  is  the  particle 
temperature  and  (u)  is  the  characteristic  energy  of  the  particle  (9:20).  This  relation 
becomes  very  important  as  many  times  the  ion  and  background  gas  temperatures  are 
assumed  to  be  equivalent.  Usually  taken  as  300°7f  or  approximately  e V  ,  Eq.  (14) 
is  used  to  evaluate  the  ion  diffusion  coefficient 
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(15) 


DP[EI  p\  =  -^Mp[EI  p\- 


Sources  and  Sinks 

Equally  important  to  the  fluid  equations  are  the  source/loss  terms.  Each  of  these 
tenns  involves  rate  parameters  that  display  the  same  energy  or  reduced  field  dependence 
seen  in  the  transport  coefficients.  These  terms  show  up  on  the  right-side  of  the  particle 
and  energy  continuity  equations  -  Eqs.  (3)  and  (5).  The  number  of  included  terms 
depends  heavily  on  the  assumptions  of  the  problem  (Appendix  A). 

The  particle  continuity  equation  may  include  such  source/loss  processes  as 
ionization,  attachment,  detachment  or  recombination.  Ionization  occurs  when  an  electron 
impacts  a  neutral  particle  above  the  ionization  threshold  energy  and  causes  that  neutral 
particle  to  lose  an  electron,  becoming  an  ion.  This  process  creates  both  an  electron  and 
an  ion  and  is  therefore  a  source  for  both  the  electron  and  ion  continuity  equations  (20:4): 

Si,e,p=neVilUe]  (16) 

In  this  equation,  vi  \ue  ]  represents  the  ionization  frequency  based  on  local  average 
electron  energy. 

Although  not  included  in  this  work,  two  processes  associated  with  an 
electronegative  gas  are  included  for  completeness.  Detachment  occurs  when  a  negative 
ion  suffers  a  collision  that  releases  its  extra  electron.  This  process  creates  one  electron 
and  one  neutral  particle  and  therefore  factors  only  into  the  electron  continuity  equation  as 
a  source: 
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(17) 


In  Eq.  (17),  kd  is  the  detachment  coefficient  and  n  represents  the  negative  ion  density 

(24:4).  This  interaction  in  reverse  is  called  attachment.  Attachment  occurs  when  an 
electron  attaches  to  a  neutral  particle  therefore  causing  the  loss  of  an  electron  and  the 
gain  of  a  negative  ion.  This  process  is  modeled  as 

^attach,-  ^ attach, e  ^e^a  0  ^) 

where  v a  represents  the  attachment  frequency  of  the  gas  (24:4). 

Finally,  recombination  covers  two  separate  processes.  The  first  occurs  when  an 
electron  recombines  with  an  ion  and  creates  a  loss  for  both  particle  species  (20:4): 

L recomb, e,p  ~  Pe-inen p  •  (19) 

where  (3 ’e_.  represents  the  electron-ion  recombination  rate.  While  this  rate  is  usually 
represented  as  a  constant  «  10  7  cm 3  Is ,  it  is  actually  dependent  on  the  local  electron 

temperature.  It  is  proportional  to  Te  72  for  lower  gas  temperatures  (ranging  from  room 

v 

temperature  up  to  the  thousands  of  Kelvin)  or  Te  72  for  higher  gas  temperatures  (9:60). 

The  second  process,  ion-ion  recombination,  occurs  when  a  negative  and  positive  ion 
collide  and  form  two  neutrals.  This  process  causes  a  loss  of  both  positive  and  negative 
ions: 


L recomb,-  ^ recomb ,p  P / 


(20) 
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Figure  11.  Ionization  Frequency  in  Helium  at  1  Torr  (17) 

Since  this  project  neglects  negative  ions,  the  only  terms  included  in  the  continuity 
source  term  will  be  ionization  and  recombination.  Combining  the  two  creates  a 
continuity  source/loss  equation  of  the  form 

Se,p  =neVi[Ue]-  Pe-inenp  (21) 

where  Se  represents  the  source/loss  term  for  both  the  ion  and  electron  continuity 

equations.  In  contrast  to  the  relatively  weak  dependence  on  mean  energy  of  the 
recombination  coefficient,  the  ionization  rate  for  Helium  gas  varies  nearly  six  orders  of 
magnitude  in  the  energy  range  5-18eV  (Fig.  11). 

The  energy  source/loss  term  S£  in  Eq.  (5)  also  has  production  and  loss  terms. 
Joule  heating  creates  a  source  (20:4): 
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(22) 


S  joule, s  —  J  '  E  —  eTe  •  E 

where  J  represents  the  electron  current  density.  Although  this  number  could  be  negative 
if  the  electron  flux  and  electric  field  are  in  the  same  direction,  it  is  more  often  a  positive 
number  and  is  therefore  normally  considered  a  source.  A  general  energy  loss  term 
creates  the  sink  (20:4): 

Lk,.,e  =neNk,[u,].  (23) 

where  kL  \ue  ]  is  the  average  energy  loss  rate  dependent  on  the  local  average  electron 
energy  (14: 1377-1378).  This  rate  represents  the  dissipation  of  electron  energy  in 
collisions  with  neutrals.  The  total  energy  source/loss  term  becomes 

Se=-ete-E-neNkL[ue].  (24) 

In  order  to  incorporate  the  energy  loss  term  into  the  energy  equation,  Eq.  (5),  the 
energy  loss  coefficient  is  expressed  in  terms  of  the  mean  electron  energy.  This  loss 
coefficient  arises  from  the  accumulated  energy  losses  due  to  elastic  and  inelastic 
collisions  with  neutrals.  The  energy  dependence  of  the  loss  coefficient  is  derived  from  a 
solution  of  the  zero  dimensional,  collisional  Boltzmann  Equation.  To  establish  the 
functional  dependence  of  this  loss  coefficient,  consider  modeling  a  homogenous,  steady- 
state  plasma  with  the  time  and  space  derivatives  in  Eq.  (5)  eliminated  and  the  balance 
expressed  as 

n,NkL[u,]=-ef,-E.  (25) 
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Assuming  the  electron  flux  to  be  only  field  dependent  and  substituting  part  1  of  Eq.  (10), 
this  balance  equation  would  be 


M,[ue\  =  -e(- /deneE)- E . 


(26) 


Rearranging  these  parameters  and  eliminating  common  terms,  Eq.  (26)  becomes 


/  \  J7 

*ik]=e(/w— • 


(27) 


Multiplying  and  dividing  the  right-hand  side  by  N  gives 


kL[ue]  =  e/ieN- 


r  EV 
KNj 


(28) 


Figure  12.  Electron  Energy  Loss  Frequency  in  Helium  (17) 
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Eq.  (28)  is  expressed  in  terms  of  electron  charge,  electron  mobility,  neutral 
number  density  and  the  reduced  field.  Since  the  mean  energy  is  only  a  function  of  the 
reduced  field,  the  functional  dependence  of  kL  can  be  inverted  to  express  this  tenn  now 
as  a  function  of  mean  energy.  The  resulting  energy  loss  frequency  expressed  as  the  loss 
coefficient  multiplied  by  the  neutral  number  density,  Nk,  \ue],  is  presented  in  Fig.  12. 

Because  the  energy  loss  rate  increases  three  orders  of  magnitude  in  the  range  of  5-18eV, 
care  must  be  taken  to  get  an  accurate  fit  or  the  numerical  results  will  be  faulty. 

Poisson’s  Equation 

While  the  fluid  equations,  transport  coefficients  and  source/loss  rates  remain  very 
important  to  the  DBD  model,  the  piece  that  ties  them  all  together  is  Poisson’s  equation. 
Ultimately,  Poisson’s  equation  relates  the  local  charge  density  to  the  electric  field  and  the 
electric  potential.  In  a  discharge,  this  relationship  governs  discharge  characteristics. 

The  relationship  starts  with  Gauss’s  Law  (21:464): 

V  •  D  =  p  (29) 

or  expanded  as: 

V  ■(sE)=q(np-ne)  (30) 

where 

E  =  -  V  (p .  (31) 
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In  these  equations,  D  represents  the  electric  displacement,  (p  is  the  electric  potential  and 
e  is  the  electric  permittivity.  It  is  important  to  note  that  the  addition  of  a  vector  symbol 
distinguishes  electric  displacement  ( D  )  from  diffusion  coefficient  ( D  ). 

It  is  common  practice  to  combine  Eqs.  (30)  and  (31)  to  generate  an  equation 
relating  the  potential  directly  to  the  particle  densities  (25:125) 

W  \sW(p)=-q{np-ne).  (32) 

Assuming  that  the  dielectric  constant  does  not  have  spatial  dependence  this  equation 
becomes 

£\/2(p  =  -q{np  ~ne).  (33) 


Eq.  (33)  is  known  as  Poisson’s  equation.  Although  Eqs.  (29),  (30)  and  (33)  are 
equivalent  this  final  transformation  ensures  that  boundary  conditions  can  be  easily 
implemented.  Electric  potential  is  a  directly  measurable  quantity  and  is  chosen  as  a 
boundary  value  both  experimentally  and  numerically.  While  some  problems  call  for  an 
electric  field  boundary  condition  (the  dielectric  boundary),  this  is  translated  into  an 
electric  potential  condition  with  Eq.  (3 1). 

With  this  definition  of  Poisson’s  equation,  each  important  feature  of  the  fluid 
approach  has  been  identified  and  developed.  The  transport  coefficients  and  source/loss 
rates  have  been  established: 

Electron  Mobility  ne  [ue  ] , 

Electron  Diffusion  De [ue]. 
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Ion  Mobility 


mp[e/p\, 


Ion  Diffusion  Dp  \E  /  p]  =  jup  [E  /  /»] , 

Electron/Ion  Source/Loss  Se  =  Sp  =  nevj  [ue ] -  /?,  ,neiip , 
Energy  Source/Loss  Se  =  -efe  ■  E  -  neNkL  [ue\. 

Each  of  the  fluid  equations  have  been  defined: 


Electron  Continuity 

Ot 

Electron  Flux 

f  =  —n  li  E  —  D  Vn  , 

e  eie  e  e  9 

Ion  Continuity 

dn  _  _ 

p  +v-rc  =sD, 

dt  p  p 

Ion  Flux 

r/;  =npJupE -DpVnp, 

Electron  Energy 

din  u  )  -  - 

v  -  e'  +  v-r=se, 

dt  £ 

Electron  Flux 

r,  =-(neue)peE-De' V{, 

Poisson’s  Equation 

V  •  (sVtp)  =  ~q{np  -  ne ). 

This  system  of  equations  can  now  be  cast  onto  a  grid  in  order  to  numerically  solve  for 
each  of  the  primary  variables  <p ,  ne,  np  and  (neue). 
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ITT-  DBD  Computational  Development 


Before  a  solution  to  this  set  of  equations  is  possible,  each  of  the  equations  must  be 
cast  into  numerical  form.  The  spatial  and  flux  discretization  schemes  as  well  as  the 
boundary  conditions  will  be  detailed.  The  time  discretization  schemes  will  be  dependent 
on  the  type  of  model  implemented. 


Table  1.  Node  and  Half  Node  Values 
Node  Values  Half-Node  Values 


(p 

E 

ne 

D 

np 

neUe 

T  T 

Ae  ’  De ,  [A  p,D  p 

£0,S 

Me,l/2’De  l/2,  A, 

Spatial  Discretization 

The  equations  are  discretized  and  cast  onto  the  staggered  grid  depicted  in  Fig. 
(13).  The  solid  circles  represent  node  locations.  At  each  node,  the  variables  appearing 
the  first  column  of  Table  1  are  determined.  The  half-node  positions,  identified  by  the  x’s 
are  shown  in  the  second  column  of  Table  1. 

Because  the  potential  and  particle  densities  will  be  used  to  specify  boundary 
conditions,  the  grid  is  set-up  so  that  there  is  a  node  for  each  boundary  and  half-nodes 
exactly  half-way  between.  In  Fig.  13,  the  shaded  nodes  represent  the  boundary  value 
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cells.  These  cells  will  not  be  included  as  part  of  the  computational  domain.  The  domain 
is  labeled  with  designators  ranging  from  0  — >  (N  - 1)  because  the  C  coding  language  uses 
this  range  for  array  indices.  As  an  example,  Fig.  14  shows  the  electric  potential  and 
electric  field  indexing  scheme  as  well  as  the  grid  spacing. 


Node 


1/2 

X 


3/2  2  (2N-5)/2 

X  •  ....  X 

1/2  Node 

Figure  13.  Grid  format 
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• 

Ei- 1/2  1/2 


-dx, 


i-1/2 


Figure  14.  Grid  Designation  -  ‘i’  Index 
In  Fig.  14,  the  grid  spacing  is  equal 

dxt  ,  =  dXj  =  dxi_in  .  (34) 

The  set-up  of  the  computational  grid  becomes  very  important  to  the  validity  of  the 
numerical  results.  The  location  and  value  of  the  boundary  conditions  is  integral  to  the 
computation.  An  accurate  distinction  of  half-cell  and  full  cell  values  in  each  of  the 
equations  is  essential  to  an  accurate  numerical  solution.  Finally,  the  choice  of 
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discretization  can  make  or  break  the  stability  of  the  numerical  technique.  One  example  is 
the  Scharfetter-Gummel  method  for  flux  discretization. 


Scharfetter-Gummel  Flux  Discretization 

Originating  in  1969  in  a  paper  addressing  Read  diodes  (26:73),  D.L.  Scharfetter 
and  H.K.  Guinmel  created  an  exponential  flux  representation  giving  greater  stability  to 
the  charged  particle  flux  calculations.  This  exponential  weighting  scheme  more 
accurately  calculates  the  flux  in  discharge  regimes  transitioning  between  field-dominated 
and  diffusion  dominated  flux.  For  DBD  applications,  the  density  flux  given  by  Eq.  (4) 
takes  the  form: 


Fs.,+i/2  =T“k ,iDs,i  exp(^s,;+i/2 )— Rv'+i^v+i  ] 

Ax  ex 


'5,1+1/ 2 


p(^+l/2  0 


(35) 


where  Ax  represents  the  grid  spacing  and 


Js,i+ 1/2 


=  -sgn[r/] 


Ms 


f+l/2 


Dsj+ 1/2 


(36) 


where  sgn[<gr]  represents  the  charge  of  the  particle  type  taking  on  a  value  of  - 1  for 
negatively  charged  particles  and  + 1  for  positively  charged  particles. 

For  greater  clarity  in  equation  development,  Eq.  (35)  will  be  redefined  with 
substitute  variables.  Taking 


HZ'J+U2. 


exp(ZM+i/2  )•  ZsM,  2 

cxp(Zi!+]/2)~l 


(37) 
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and 


F2[Z,MI  2]  = 


ex 


7 

Z"i,i+l/2 

P(^j,/+l/2 


(38) 


then  Eq.  (35)  becomes 


r,,m/2  —  A  ^s,!+i^s,!+i^,2[ZijJ.+1/2]). 


Ax 


(39) 


Applying  this  same  procedure  to  the  energy  flux,  Eq.  (6)  takes  the  form 


r,  =  7-  ((«,«,),  De,iFl[z  1/2  ]  -  );+l  De .  <+1F2[z, e><+1/2  ]) . 

Ax 


(40) 


Eqs.  (39)  and  (40)  will  be  used  for  all  fluxes  in  the  domain  except  the  boundary  flux. 


z 


Figure  15.  Fl[z]  and  Fl[z] 
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The  utility  of  the  Scharfetter-Gummel  flux  discretization  method  can  be  better 
described  through  example.  Fig.  15  shows  that  Fl[z]  and  F2[z\  are  symmetric  across 
the  y-axis.  Manipulating  the  discretized  form  of  Eq.  (31)  yields 

(<Pm  -<Pi)  =  -E-Ax.  (41) 

Replacing  the  electric  potential  tenn  in  Eq.  (36)  we  can  relate  Z  directly  to  E  as 


7  —  „  Fs,i+ 1/2  p 

Js,i+l/2  ~S~^  ^i+l/2 


■  Ax . 


(42) 


In  the  diffusion  limit  of  the  flux  approximation,  the  field  contribution  is 
essentially  zero  (27:9).  From  Eq.  (42)  this  means  that  Z  will  be  approximately  zero.  In 
this  limit,  the  F-terms  become 

Z,.„2  *0->(Fl[zw,,]»F2[zw,2]»l).  (43) 


Now  the  flux  equation  becomes 


ri+i/2  =7 ~(niDi  ~nMDM). 
Ax 


(44) 


Assuming  a  spatially  constant  diffusion  coefficient,  it  can  be  cast  into  the  familiar 
equation  seen  in  part  2  of  Eq.  (10): 


r,+i/2  —  D 


Ax 


(45) 
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The  drift  limit  of  the  Scharfetter-Gummel  flux  can  be  detennined  in  a  similar 
manner.  In  this  limit,  the  electric  field  would  be  either  a  large  negative  or  large  positive 
value.  The  value  of  Z  also  becomes  either  significantly  greater  than  or  less  than  1 .  In 
this  regime,  one  of  the  F-terms  will  go  to  zero  (27:9): 


D 


sgn[g]-fi  »  — 1 


+1/2 


Mi 


■^1+1/2  >>  1 


+1/2 


>1  [z«,d->z, 


+1/2 


(46) 


or 


sgn[q\-E  «  D‘+V2  ->  Z.  , ,,  «  1 


A 


1+1/2 


1+1/2 


Fl[Zi+ 1/2]^  0  ' 

F’2[Z(.+1/2]  — >  |Z<+1/2|, 


(47) 


Combining  Eqs.  (39),  (42)  and  (46),  the  flux  becomes 


^•+1/2  =  sgn[^]- 

niMi+\/2^i 


+1/2  * 


(48) 


Combining  Eqs.  (39),  (42)  and  (47),  the  flux  would  be 


r/+i/2  =-sgn[g]- «,.+1//, +1/2^, 


+1/2  • 


(49) 


Eqs.  (48)  and  (49)  also  help  relate  how  the  Scharfetter-Gummel  flux  discretization 
is  a  fonn  of  upwinding  (27:10).  The  conditions  of  Eq.  (48)  prescribe  that 
sgn[g]  •  Em/2  »  DMI2  / ni+i/2  ■  This  means  that  the  field  must  be  large  and  the  charge  of 

the  particle  and  the  polarity  of  the  field  are  the  same.  This  applies  to  a  negative  electric 
field  for  electrons  and  a  positive  electric  field  for  ions.  Under  these  conditions,  the  flux 
at  the  /  +1/2  half-node  is  based  on  the  particle  density  at  node  i .  If,  on  the  other  hand, 


32 


the  charge  of  the  particle  and  polarity  of  the  field  are  opposite,  the  flux  at  the  i  + 1/2 
half-node  is  based  on  the  particle  density  at  node  /  +  1 .  The  upwinding  process 
contributes  to  the  stability  of  the  Scharfetter-Gummel  flux  discretization  method.  This 
method  is  used  to  calculate  all  fluxes  except  the  boundary  values. 

Boundary  Conditions 

For  this  project,  the  flux,  potential  and  particle  densities  are  specified  at  the 
boundaries.  The  boundary  particle  densities  will  be  treated  in  two  different  ways 
depending  on  the  type  of  model  being  implemented.  In  some  cases,  the  boundary 
densities  were  taken  to  be  zero  because  the  boundary  cells  do  not  factor  into  the 
continuity  or  flux  calculations.  In  other  cases,  it  is  necessary  to  set  the  ion  density 
gradient  equal  to  zero  at  the  boundary  so  that  the  ion  boundary  flux  is  field  driven  only 
(14: 1379).  The  particular  cases  to  which  each  specification  applies  will  be  covered  in  the 
next  chapter. 

For  the  electron  and  energy  flux  at  the  boundaries,  this  code  implements  a  thermal 
flux  with  secondary  emission  for  all  exposed  electrodes.  For  the  electrons,  the  boundary 
thermal  flux  becomes 


^e,H2  — ,N-1-1I2  I-  ^  neve,th 


where  ve  th  is  the  electron  thermal  velocity  which  can  be  calculated  as 


v 


e,th 


n  ■  me 


(50) 


(51) 
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me  is  the  electron  mass  and  T,  is  the  electron  temperature  (14:1379).  The  thermal  flux  is 
always  directed  toward  the  electrode  so  that  the  sign  in  Eq.  (50)  will  be  negative  for  the 
left  electrode  and  positive  for  the  right.  Substituting  known  values  and  relating  kBTe 
back  to  the  electron  average  energy  this  equation  becomes 

ve,th  ~ 4.19x  105^/m7 .  (52) 


The  electron  energy  flux  takes  a  very  similar  form  (14:1379): 


r.4,2  —^e,N-\-V 2  +"  ^  neVe,th(jiB^1e)  ■ 


(53) 


Again  substituting  values,  this  equation  relates  to  the  electron  energy  density  as 


r  =r 

A  £,1/2  x  e,N- 1-1/2 


=  +\{neUe)VeJh- 


(54) 


The  flux  from  secondary  emission  captures  the  flux  of  electrons  being  ejected 
back  into  the  plasma  due  to  ion-electrode  impact.  This  flux  is  characterized  by 


Te  =  -JT, 


(55) 


where  y  is  defined  as  the  secondary  emission  coefficient.  This  coefficient  is  usually 
modeled  as  a  constant  with  values  ranging  between  0.0  and  0.5  (16:6). 

The  ion  flux  at  the  boundaries  is  modeled  as  the  Scharfetter-Gummel 
representation  of  the  flux  at  that  cell.  The  ion  flux  to  the  electrodes  is  field-driven  when 


34 


the  drift  velocity  is  directed  toward  the  wall  and  zero  otherwise  (14: 1379).  The  ion  flux 
toward  the  electrode  becomes 


r„,,„  =  T(F1[Z|;j  ]_  F2[z„3  ])Bp  , 


(56) 


to  the  left  electrode  and 


^p,N-i/2  ~  ^  (^[Xv-3/2  ]  ^2[Z  v_3/2  J)n p^N_2D p 


N- 2 


(57) 


to  the  right. 

For  the  potential,  the  boundary  conditions  simply  specify  the  value  of  the 
potential  at  the  electrode.  These  values  vary  greatly  depending  on  which  test  case  is 
being  implemented.  For  each  case,  however,  the  potential  at  the  electrodes  is  held 
constant  for  each  time  step.  If  a  dielectric  is  covering  the  electrode,  Poisson’s  calculation 
must  take  into  account  the  charge  build-up  on  the  dielectric  surface. 

The  Dielectric  Boundary 

The  dielectric  configuration  used  for  this  project  is  shown  in  Fig.  16.  Both 
electrodes  and  the  dielectric  barrier  are  considered  to  be  infinite  in  transverse  extent. 
While  there  are  several  models  that  could  be  used  to  address  the  charge  build-up  on  the 
dielectric  surface  (6:313;  28: 167;  29:98),  this  project  will  use  the  method  prescribed  by 
Boeuf  (14: 1379).  The  dielectric  surface  is  modeled  as  sticky  accumulating  all  charge 
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striking  its  surface.  Additionally,  this  model  assumes  that  once  the  charges  strike  the 
surface,  any  ion  and  electron  pair  recombine  instantaneously. 


t 


Electrode 


1 


Dielectric.  \ 


Electric  field 


Infinite  extent 


1 


Electrode 


Figure  16.  1-D  Dielectric  Configuration 

The  geometry  of  the  dielectric  and  impinging  flux  is  represented  by  Fig.  17.  The 
surface  charge  accumulates  as 


ds 


8t 


=  elr 


(58) 


where  crds  is  the  surface  charge  density  accumulated  on  the  dielectric.  Because  the 
charges  stick  to  the  dielectric  surface,  they  only  contribute  to  the  surface  charge  density 
and  not  the  individual  particle  densities  ne  and  np  .  The  particle  nds  are  zero  and  the 

electron  flux  onto  the  surface  is  taken  as  the  thermal  flux  without  secondary  emission 


e,ds- 1/2 


le,ds- 1  v  e,ds—\/  2,th  * 


(59) 


The  ion  flux  is 
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^p,ds— 1/2  ^ ^  ])  * 


(60) 


The  charge  build-up  is  ignored  in  the  flux  equations  but  accounted  for  in 
Poisson’s.  At  the  dielectric  surface,  Poisson’s  equation  is  altered  to  account  for  the 
surface  charge.  In  Fig.  17,  the  permittivity  of  the  dielectric  substance  sx  modeled  to  the 
right  of  (pds  and  left  of  cpds+l  and  the  permittivity  of  free  space  s0  is  modeled  to  the  left 
of  (pds .  Starting  with  the  discretized  form  of  Eq.  (30) 


sxEt 


ds+ 1/2 


Eds- 


Ax 


ds 


Ax 


—  =  q{np,ds  ~ne,ds)- 


ds- 1 


(61) 


Upon  substituting  the  potentials  and  rearranging 


(gl^Ufc-l  +  ZpAXds  +  ZpAXdstPds-l 

^ds^ds-l 


=  —q{n 


p,ds  ^ e,ds  , 


(62) 


gives  Poisson’s  equation  accounting  for  the  surface  charge  density  on  the  dielectric 
surface.  Now  the  electric  potential  can  be  found  throughout  the  computational  domain, 
to  include  the  dielectric  cells,  by  one  of  the  methods  addressed  in  the  Numerical  Methods 
section. 


X 


X 


ds+1/2 


•  X  • 

ds+1 


I - dxds- 1 - 1 - dXds - 1 

Figure  17.  Dielectric  Grid  Geometry 
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1-D  Numerical  Methods 


Two  different  numerical  models  are  used  for  the  calculations  in  1-D.  The  first 
model  involves  the  sequential  solution  to  Poisson’s  equation,  the  electron  and  ion 
continuity  and  the  electron  energy  continuity  equations  -  Eqs.  (3),  (5)  and  (33) 
respectively.  The  second  involves  the  simultaneous  solution  to  all  of  these  equations. 
Each  of  these  methods  will  be  covered  in  detail. 

Semi-implicit  Sequential 

The  first  numerical  model  involves  a  semi-implicit  sequential  solution  to 
Poisson’s  equation  and  the  electron,  ion  and  energy  continuity  equations.  The  semi- 
implicit  designation  indicates  that  some  variables  in  each  equation  will  be  evaluated  at 
the  previous  time  step  and  others  at  the  current  time  step.  In  finite  difference  form  using 
‘k’  as  the  time  index  and  ‘i’  as  the  spatial  index,  Poisson’s  equation  (non-changing 
dielectric)  becomes 


jpiA  -2(p,  +(pM  j  k  k\ 

£ — - — ^ —  -  -v\nP.i  ~  "«.<;• 


(Av)2 


(63) 


Eq.  (62)  with  a  time  index  ‘k’  on  all  potential  and  density  variables  is  the  form  for  the 
dielectric  cell.  The  electron  and  ion  continuity  equations  become 


k+ 1  k  t"»&+1  t~'&+1 

^ e,i  ^ e,i  ^e,(+l/2  ^  e,i— 1/2  k+\  k 

-H - n„  ■  v  . 


At 


Ax 


k*,]+ =° 


(64) 


and 
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(65) 


lk+l 

p,l 


—  n 


p,i 


At 


+  ■ 


^k+ 1 
/?,/+!/  2 


-r 


t+i 

p,i-l/2 


Ax 


+  fin 


k+ 1  A+l 
”p./ 


=/i*rv* 


The  source  term  in  Eq.  (65)  must  be  on  the  right-hand  side  because  it  involves  the 
electron  density.  This  term  is  still  at  the  advanced  time  step,  however,  because  the 
electron  continuity  calculation  is  perfonned  first.  The  flux  representation  for  Eqs.  (64) 
and  (65)  is 


T"*  k+ 1 

1  e(,p),i+ 1/2 


Ax 


(n 


k+ 1 
e(,p),i 


k 

e(,p),i+l 


/  2 


k+\ 

e(,p),i+ 1 


Z) 


k 

e(,p\i+\ 


k 

e(,p),i+ 1/2 


1) 


(66) 


where 


zg(.p),-+i/2  =-sgnb]^;('/,)-,+1/3  (?4  -?>*). 


D 


(67) 


<?(,/?),/+!/ 2 


Lastly,  the  energy  density  continuity  equation  looks  like 


\neUe)i  —\neUe)i  ,  ^  s.i+1/2  _  ^  £,(-1/2  f=.A+l  £ ,A  A+l  Ar,  [  A  1 

- S - + - At - =  “r"  E‘  J 


(68) 


where 


t-<  k+\ 

1  e,i+ 1/2 


Ax 


+1/2 


(69) 


Note  that  the  assumption  that  all  energy  quantities  are  given  in  eV  eliminates  the 
electron  charge  e  originally  in  Eq.  (68). 
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Figure  18.  Semi-Implicit  Sequential  Iterative  Scheme 

The  sequential  numerical  method  solves  Eqs.  (63),  (64),  (65)  and  (68)  in  succession. 
These  solutions  give  the  updated  values  for  (p,ne,  np  and  (, neue )  respectively.  Fig.  18 

represents  the  iterative  scheme  where  one  iteration  involves  the  consecutive  solutions  to 
Poisson’s,  electron  continuity,  ion  continuity  and  electron  energy  equations. 

The  iterative  scheme  for  the  sequential  method  is  shown  in  Fig.  18.  One  iteration 
involves  the  consecutive  solutions  Eqs.  (63),  (64),  (65)  and  (68).  Poisson’s  equation  is 
solved  using  an  SOR  routine  while  each  of  the  continuity  equations  is  solved  using  a 
tridiagonal  solver.  After  each  iteration,  the  reduced  field/energy-based  rates  are  reset. 
The  electron  mobility  and  diffusion  coefficients  as  well  as  the  ionization  and  energy  loss 
rates  are  set  using  the  updated  average  electron  energy.  The  ion  mobility  coefficient  is 
adjusted  with  the  updated  reduced  field  E  /  p  and  the  ion  diffusion  coefficient  is  set 
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using  Eq.  (15)  assuming  the  ion  temperature  is  approximately  equal  to  the  gas 
temperature  at  0.02eV  .  Once  these  parameters  are  updated,  the  code  recalculates  the 
time  step  using  Eq.  (70)  and  returns  to  solving  the  primary  equations. 

For  each  iteration,  one  time  step  is  taken  where  each  time  step  is  limited  by  the 
dielectric  relaxation  time  tdi  where 


tdi  = 


e\neMe+nUL 


7 


(70) 


This  limitation  factors  in  because  the  sequential  solution  can  not  overstep  the  movement 
of  each  of  the  particles  in  the  continuity  equation.  The  drift  portion  of  the  flux  equations 
describes  the  motion  of  charged  particles  in  an  electric  field.  Yet  it  doesn’t  take  into 
account  variations  in  the  electric  field  as  each  of  the  particles  moves  (15:5609).  If  a  time 
step  is  taken  that  exceeds  the  dielectric  relaxation  time,  the  particles  move  enough  to 
significantly  alter  the  electric  field  and  the  field  driving  their  motion  becomes  invalid.  If 
this  occurs  too  often,  the  code  quickly  becomes  unstable  causing  either  a  crash  or 
inaccurate  results.  The  key  to  stability  lies  in  an  accurate,  updated  solution  to  Poisson’s 
equation  using  the  SOR  algorithm. 

SOR 

The  method  of  successive  overrelaxation  (SOR)  is  based  on  the  Gauss-Seidel 
iterative  algorithm.  Both  algorithms  are  explained  in  depth  in  several  computational  texts 
(30:863-869;  31:192-194;  32:162-166).  For  this  project,  the  SOR  algorithm  was 
implemented  as  a  root  solving  routine  where  Eq.  (63)  takes  the  form  of 
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(71) 


2  <P*  +  <Pi+ 1  +  ~ 
— - '  £ 


Assuming  that  we  currently  know  all  of  the  variables  in  part  2  of  Eq.  (71),  the  goal  is  to 
find  out  what  values  of  (pk  , ,  cpk  ,  and  (p'l ,  make  this  equation  true  for  all  spatial  grid 
points  i  =  1,2, . . . ,  N  -  2  .  The  first  step  is  to  guess  values  for  the  potential  for  the  first 
iteration,  call  these  values  (p°,d'k .  Next,  solve  for  an  updated  (p*,k  using  these  guess 
values: 


„ old,k 
<Pi-\ 


<i-q- 

£ 


(Ax)2  (nkp  i  -  nke  i ) 


(72) 


Now  take  a  weighted  sum  of  the  old  and  the  updated  potentials  and  designate  the  new 
value  of  the  electric  potential: 

(Pi  ~  0)(P,  +  (1  -  co)(pl  (73) 


where  co  is  the  weighting  parameter.  The  value  of  the  weighting  parameter  that  gives  the 
fastest  convergence  (32:166-167)  is  found  by  solving 


(74) 
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where  coopl  is  the  optimal  weighting  parameter  to  achieve  the  fastest  convergence. 
Finally,  check  the  convergence.  For  this  code,  the  convergence  criteria  was  g  <  le  -  5 
where 


Cycle  through  Eqs.  (72),  (73)  and  (75)  for  all  spatial  values  until  the  convergence  criteria 
is  reached.  At  this  point,  the  electric  potential  contained  in  the  ‘new’  array  is  taken  to  be 
the  value  of  the  potential  for  the  current  time  step.  This  potential  is  used  to  solve  for  the 
electric  field  which  is  used  for  the  remainder  of  the  semi-implicit  iteration  to  solve  the 
fluid  equations.  For  the  semi-implicit  method,  all  the  potentials  are  values  from  the 
previous  iteration  because  they  are  simply  accounting  for  the  density  changes  that 
occurred  during  the  last  cycle. 

Tridiagonal 

The  fluid  equations  will  be  solved  using  a  generalized  Thomas  Algorithm  (33). 
This  algorithm  is  commonly  used  as  an  efficient  solver  for  systems  of  equations  that  can 
be  cast  into  tridiagonal  form.  With  a  little  manipulation,  all  of  the  fluid  equations  can  be 
formatted  as  tridiagonal  systems  of  equations.  Using  the  electrons  as  an  example,  first 
combine  Eqs.  (64)  and  (66)  to  give 
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nk+l  —nk- 

e,i  e,i 


At 


+ 


Ax 


—  (nk+\Dk.  ,F\\zk.  ,n]-nk+lDk.Fl\zk.  ,,,1) 

*  \  e,i— 1  e,i— 1  L  e,i-\/2  I  e,i  e.i  L  e,/— 1/2  1/ 

Ax 


•  (76) 


Ax 


*+l  IT  k  1  .  /o  &+1  &  rv 
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Multiplying  through  by  At  ■  (Ax)  and  rearranging  terms,  this  equation  becomes 


At  •D«*,_,Fl[zi_„JJ)+ 

A (Ar)2  +A(D,‘,Fl[z,*,.l(2]+Ar.£l,‘,F2fc;H,2l 

“  -('Vc)'  A/  ■  v'.': [«'  ,]+  ('\,)  A/  ■  fln^  . 


+ 


(Ax)2 


(77) 


This  equation  is  now  in  tridiagonal  form.  Taking  the  coefficients  of  the  discretized 
electron  densities  at  the  node  locations  ne  M ,  ne .  and  ne  M  and  associating  them  with  / , 

dj  and  a  .  ,  then  designating  the  known  variables  of  the  right-hand  side  to  equal  /'  this 
equation  looks  like 


<»  )+ C.  (“-)=/.  • 


(78) 


The  coupled  equations  can  take  on  the  matrix  form  Ax  =  b: 
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(79) 
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In  this  equation,  A  is  the  (N  -  2)x  (N  -  2)  coefficient  matrix  for  the  electron  density  at 
the  current  time  step,  x  is  the  vector  of  unknown  electron  densities  at  the  current  time 
step  and  b  is  vector  of  known  values  consisting  of  the  electron  density  at  the  previous 
time  step  and  the  ionization  source  term.  Now  the  generalized  Thomas  algorithm  can  be 
easily  implemented  to  solve  this  system  of  equations  for  the  electron  density  at  the  next 
time  step. 

The  ion  and  energy  continuity  equations  are  manipulated  in  the  same  manner  used 
to  create  Eq.  (77).  Moving  the  source  term  to  the  right-hand  side  and  changing  all  other 
particle  densities  to  ions,  the  ion  continuity  equation  is  identical  to  Eq.  (77).  The  energy 
density  becomes 
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Using  the  new  values  for  the  lt ,  dj ,  u ,  and  f.  coefficients,  the  same  tridiagonal  method 

is  used  to  solve  for  the  advanced  time  step  ion  and  energy  densities. 

Semi-Implicit  Summary 

Overall,  the  semi-implicit  sequential  method  is  relatively  easy  to  code  and  gives 
accurate  results  as  long  as  the  time  step  remains  within  the  dielectric  relaxation  time. 

This  time  step,  however,  was  detennined  to  be  a  constraint  for  goals  of  this  project.  As 
Eq.  (70)  relates,  the  dielectric  relaxation  time  is  inversely  proportional  to  the  particle 
densities.  This  time  step  significantly  diminishes  with  nearly  any  increase  in  the  charge 
densities. 

As  an  example,  helium  gas  at  1  torr  with  an  assumed  initial  charge  density  of 
1x10 15  m  3  (an  ionization  fraction  on  the  order  of  1  x  10  7 )  has  a  dielectric  relaxation  time 
«  0. 1 5  ns  before  any  growth  due  to  ionization  is  taken  into  account.  If  the  goal  is  to 
cover  one  cycle  of  an  RF  discharge  with  a  10MHz  frequency,  it  would  only  require 
around  670  iterations.  If  the  goal  is  a  steady-state  solution  (-500  cycles),  however,  this 
time  step  becomes  a  serious  hindrance.  This  is  why  the  fully-implicit  simultaneous 
numerical  method  is  the  next  step.  This  method  allows  a  larger  time  step  which  greatly 
enhances  the  modeling  power  of  the  code. 

Fully  Implicit  Simultaneous 

The  second  numerical  model  involves  a  fully  implicit,  simultaneous  solution  to 
Poisson’s  equation  and  the  electron,  ion  and  energy  continuity  equations.  ‘Fully  implicit’ 
indicates  that  the  values  for  all  potentials,  particle  densities,  energy  densities  and 
source/loss  rates  will  be  evaluated  at  the  advanced  time.  This  algorithm  becomes  more 
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involved  than  the  sequential  method  because  it  must  address  the  non-linear  nature  of  the 
coupled  system  of  equations. 

The  finite  difference  fonn  of  these  equations  utilizes  the  Crank-Nicholson  method 
for  stabilization  and  the  Newton-Raphson  method  to  linearize  and  solve  the  system.  For 
this  section,  the  designation  of  the  iteration  and  time  variables  are  changed  significantly. 
An  ‘in’  designates  the  previous  time  step,  ‘m+1’  the  current  time  step,  ‘k’  the  previous 
iteration  and  ‘k+1’  the  current  iteration. 

Crank-Nicholson 

The  Crank-Nicholson  method  is  an  implicit  technique  that  provides  second  order 
accuracy  in  both  time  and  space  dimensions  (34:841).  This  method  uses  the  average  of 
the  flux  at  the  previous  time  step  m  and  the  current  time  step  m  + 1  to  replace  the 
divergence  tenn  in  the  continuity  equations.  The  Crank-Nicholson  differenced  forms  of 
Eqs.  (3)  and  (5)  become 
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As  the  time  step  is  increased,  this  averaging  method  adds  needed  stability  to  the  system 
of  highly  non-linear  equations.  In  order  to  solve  the  equations,  the  Newton-Raphson 
algorithm  is  employed. 

Newton-Raphson 

In  its  simplest  form,  the  Newton-Raphson  method  is  a  root-finding  algorithm 
(34:147-151).  For  the  purpose  of  the  project,  however,  it  is  used  as  a  technique  to  solve 
systems  of  non-linear  equations.  While  McGrath  (35)  gives  an  excellent  introduction  to 
the  Newton-Raphson  method  as  well  as  a  simple,  clear  example  of  how  it  can  be  used, 
some  of  the  basics  will  be  covered  as  they  apply  to  this  project. 

Before  beginning  a  discussion,  it  is  important  to  recall  some  of  the  basic  attributes 
of  the  system  of  equations.  The  four  primary  equations  of  interest  are  Poisson’s  equation, 
the  electron  continuity,  ion  continuity  and  energy  continuity  equations.  The  four 
variables  are  (pn  nei ,  npi  and  {neue)r 

The  Newton-Raphson  algorithm  starts  by  assuming  that  the  variable  at  the  next 
time  step  will  be  equal  to  the  value  of  the  variable  at  the  current  iteration  plus  some  delta 
value.  Using  the  four  variables  as  examples,  this  means  that 
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(83) 


where  the  5  terms  are  the  corrections  to  the  variables  for  this  time  step  but  the  previous 
iteration.  As  a  quick  reminder,  k  + 1,  m  + 1  indicates  the  current  time  step  and  current 
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iteration,  k,  m  + 1  indicates  the  current  time  step  and  previous  iteration.  Any  variable  at 


k  +  1,  m  +  1  can  be  approximated  in  this  way  as  long  as  the  variable  is  linear.  The  process 
becomes  slightly  more  complicated  for  the  non-linear  variables. 

Using  the  electron  flux  at  the  ‘i+1/2’  spatial  location  as  an  example,  it  is  easy  to 
see  from  Eqs.  (36)  and  (39)  that  the  flux  is  dependent  on  the  electron  density  and  the 
electric  potential  at  ‘i’  and  ‘i+1’  spatial  locations.  For  this  equation,  the  changes  to  the 
flux  must  account  for  these  four  dependencies.  The  flux  correction  is  still  represented  as 
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In  this  equation  the,  the  total  8  for  the  flux  is  the  sum  of  the  8  ’s  of  each  of  the 
constituent  variables  times  their  partial  derivatives.  This  same  method  is  used  for  all  the 
non-linear  terms  in  the  equations. 

Now  using  the  electron  continuity  equation  as  an  example,  the  entire  Newton- 
Raphson  development  will  be  covered.  To  begin,  Eq.  (3)  will  be  written  in  implicit  form 
with  a  Crank-Nicholson  representation  of  the  flux  tenn: 
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Moving  each  of  the  variables  to  the  left-hand  side,  this  equation  can  be  cast  into  root 
solving  fonn.  Multiplying  by  At ,  Eq.  (86)  becomes 
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where  NE[<p,  ne,  nesi  \  represents  the  total  electron  equation  in  root-solving  form.  Next, 

all  terms  at  the  current  time  step  and  current  iteration  are  substituted  with  their  linear 
equivalents: 
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Understanding  that  both  flux  tenns  and  the  ionization  rate  are  non-linear  tenns,  in 
addition  to  Eq.  (85)  the  following  values  will  be  substituted: 
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and 
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Simplifying  and  rearranging,  this  equation  becomes 
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This  process  is  repeated  for  the  remaining  three  equations.  The  original  and  final 
Newton-Raphson  form  of  the  remaining  three  equations  can  be  found  in  Appendix  B. 
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The  system  of  equations  is  now  linear  and  can  be  cast  into  the  matrix  form 


Ax  =  b  .  This  matrix  equation  looks  like: 
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where  An ,  Bn  and  Cn  are  4x4  sub-matrices,  or  blocks,  consisting  of  the  coefficients  of 
the  S  terms,  xn  is  a  four-tenn  sub-vector  of  the  unknown  5  ’s  and  bn  is  a  four-term  sub¬ 
vector  of  the  variables  not  involving  a  5  .  Looking  at  just  one  submatrix  equation  shows 
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where  PHI ,  NE ,  NP  and  EN  represent  the  total  root-solving  representation  of  each  of 
the  variable  equations  and  PHl(s),  Ne(s),  Np{§ )  and  En[s^  represent  the  terms  in 
each  equation  that  do  not  involve  delta  values.  For  each  of  the  coefficients  of  the  sub¬ 
matrices  as  well  as  the  non-delta  terms,  see  Appendix  B.  Now  the  Newton-Raphson 
linearization  is  complete. 
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From  here,  the  problem  becomes  an  iterative  root  solving  routine  for  a  system  of 
linear  equations.  For  each  time  step,  the  total  coefficient  matrix  is  solved  for  the 
unknown  8 ’s.  These  8  ’s  are  added  to  the  variables  at  the  current  time  step  but  previous 
iteration.  This  updates  the  variables  for  the  current  time  step  and  current  iteration: 
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Ideally,  the  roots  of  the  primary  equations  are  found  when  all  of  the  8  ’s  go  to 
zero.  It  is  considered  adequate,  however,  that  all  of  the  8  values  be  less  than  some 
tolerance  that  is  significantly  less  than  one.  For  this  project  taking  each  of  the  variables 
as  x ,  the  tolerance  was  g  <  \e  -  5  where 
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When  this  condition  is  met,  the  next  time  step  is  taken.  The  last  converged  solution  is 
used  as  an  initial  guess  for  the  next  solution.  This  root  solving  routine  is  repeated  until 
the  desired  total  time  or  number  of  time  iterations  is  reached. 

Fully  Implicit  Summary 

When  all  of  the  primary  equations  are  solved  simultaneously,  the  time  step 
constraint  given  by  Eq.  (70)  relaxes  significantly.  Each  of  the  variables  takes  the  same 
time  step  concurrently  and  there  is  no  particle  motion  that  is  not  accounted  for  by 
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Poisson’s  equation.  Additionally,  the  Crank-Nicholson  method  adds  stability  and  allows 
for  an  even  greater  time  step.  For  the  ambipolar  test  case,  for  example,  while  the  semi- 
implicit  method  was  limited  to  time  steps  between  0.1  -  Ins  ,  the  fully  implicit  method 
allowed  a  time  step  on  the  order  of  0.  \/js  ,  two  to  three  orders  of  magnitude  greater. 

Still  this  algorithm  is  not  without  limitations.  While  the  Newton-Raphson  method 
provides  a  powerful  tool  for  linearization  and  solution  of  these  systems  of  equations,  it  is 
limited  by  the  requirement  that  the  initial  guess  be  adequately  close  to  the  true  solution. 
Typically,  the  solution  on  the  previous  time  step  is  used  as  an  initial  guess  for  the  current 
time  step.  In  order  for  the  solution  to  converge,  there  still  must  be  some  limit  on  the  time 
step  based  on  the  needed  accuracy  of  the  initial  guess. 

This  is  especially  true  for  modeling  cases  where  the  sheath  region  plays  a 
significant  role.  In  these  test  cases,  the  electric  field  is  large  near  the  cathode  and  the 
electron  particle  densities  change  quickly.  If  the  time  step  is  too  large,  the  conditions 
change  too  much  for  the  previous  solution  to  provide  a  satisfactory  guess  for  the  next 
solution.  While  there  is  still  an  increase  in  the  allowed  time  step,  it  is  not  as  significant  as 
the  ambipolar  case.  Each  of  these  test  cases  is  covered  in  detail  in  the  following  chapter. 
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IV.  Computational  Validation  and  DBD  Results 


In  order  to  accurately  characterize  the  DBD  in  one-dimension,  the  numerical 
model  must  be  extensively  tested  to  ensure  that  the  code  produces  expected  results.  The 
semi-implicit  numerical  model  was  checked  against  a  variety  of  analytic  and  previously 
validated  test  cases  before  the  DBD  characterization  tests  were  run.  Because  of 
difficulties  implementing  the  energy  equation  into  the  fully  implicit  model,  only  the 
ambipolar  test  case  was  accomplished.  In  each  test,  the  spatial  domain  was  divided  into 
101  cells  -  99  computational  and  2  electrode  boundary  cells.  Unless  otherwise  specified, 
the  length  of  the  cavity  is  0.04m,  the  electron  boundary  flux  is  thermal  with  secondary 
emission  according  to  Eqs.  (52)  and  (55)  and  energy  density  boundary  flux  is  thermal 
according  to  Eq.  (54). 


Validation 

The  validation  test  cases  include  an  analytic  examination  of  Poisson’s  equation, 
transient  sheath  analysis,  ambipolar  analytic  comparison  as  well  as  a  radio-frequency 
comparison  to  previously  validated  results.  Each  of  these  tests  is  covered  in  detail  below. 
Several  of  the  analytic  comparisons  are  listed  in  terms  of  relative  error.  Taking  ‘x’  to  be 
the  variable  of  interest  this  relative  error  would  be: 


%error  = 


x 


analytic 


-x 


computational 


X 


analytic 


x  100% 


(96) 
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If  the  results  are  in  array  fonnat  and  will  be  compared  across  an  entire  domain,  the  errors 


will  be  the  average  %  error  where 


avg  _  %error 


N-2 

I 


X  —  X 

analytic  computational 


v  analytic 


x  100% 


N-2 


(97) 


Analytic  Poisson 

An  accurate  solution  to  the  DBD  problem  hinges  on  an  accurate  solution  to 
Poisson’s  equation.  The  first  test  case  was  a  comparison  of  the  computational  results 
with  an  analytic  solution  of  Poisson’s  equation.  This  test  was  run  for  both  the  semi- 
implicit  and  fully  implicit  numerical  schemes  and  included  only  the  Poisson  solver. 

First,  all  charge  densities  were  eliminated  and  a  potential  was  applied  to  the  right 
and  left  electrodes.  For  this  test  case,  the  analytic  solution  shows  a  linear  dependence  of 
the  electric  potential  on  x: 


<p(x) =<Pl  ^  x  +  P0  (98) 

where  (p0  is  the  electric  potential  at  the  left  electrode,  <pL  is  and  electric  potential  at  the 

right  electrode,  L  is  the  electrode  separation  and  x  is  the  domain  location  referenced 
from  the  left  electrode  for  this  project.  Test  cases  were  run  with  -  500F  ,  -  10F  and 
100F  applied  to  the  left  electrode  and  respectively  OF,  10F  and  -305F  applied  to  the 
right  electrode.  Both  methods  quickly  reached  the  convergence  criteria  of  1 .0e  -  5  . 
Convergence  was  determined  by  Eqs.  (72)  and  (95)  for  the  semi-implicit  and  fully 
implicit  methods  respectively. 
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Next,  a  uniform  charge  density  is  distributed  throughout  the  computational 
domain.  Assuming  this  charge  density  is  some  px ,  the  analytic  solution  to  this  test  case 
is  the  simple  quadratic 


pM=a*2+— — V^^o-  (") 

The  same  three  voltage  configurations  were  tested  with  net  positive  and  negative  charge 
densities.  Again,  both  solvers  quickly  reached  the  convergence  criteria.  The  solution  to 
the  first  test  case  for  both  the  zero  and  constant  charge  densities  are  shown  in  Fig.  19. 

The  final  check  on  the  Poisson  solvers  involved  the  insertion  of  various  dielectrics 
each  with  a  different  permittivity.  Once  again,  a  test  with  no  charge  density  and  constant 
charge  density  was  performed.  These  test  cases  were  checking  for  the  continuity  of  D  at 
the  dielectric  surface: 


Dx=D2  (100) 

or 

£xEx=£2E2  (101) 

where  the  ‘1’  subscript  will  designate  quantities  to  the  left  and  a  subscript  ‘2’  will 
designate  quantities  to  the  right  of  the  dielectric  change.  Both  the  semi-implicit  and  fully 
implicit  solvers  again  quickly  reached  the  convergence  criteria  and  showed  an  exactly 
continuous  electric  displacements  at  the  boundary.  Fig.  20  shows  the  results  of  both  the 
zero  charge  density  and  constant  charge  density  test  cases.  The  test  case  shown  involves 
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Electric  Potential  (V)  Electric  Potential  (V) 


electric  potentials  of  -  500V  applied  to  the  left  electrode  and  OF  applied  to  the  right 


electrode,  a  pennittivity  of  2 e0  to  the  left  and  a  permittivity  of  0.5^0  to  the  right. 


Figure  19.  Poisson  Check  -  No  Dielectric 
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With  these  tests  completed  both  the  fully  implicit  and  semi-implicit  Poisson 
solvers  with  and  without  the  dielectric  have  been  validated.  This  integral  piece  of  the 
code  can  now  be  used  in  the  remainder  of  the  test  cases  without  concern  for  its  accuracy. 
This  accuracy  becomes  extremely  important  for  the  transient  sheath  test  case  where  the 
electric  potential  changes  very  quickly  in  the  sheath  region. 

Transient  Sheath 

The  transient  sheath  problem  served  as  an  excellent  initial  test  of  the  semi-implicit 
solver  two  reasons.  First,  it  creates  a  region  with  a  strong  electric  field  (sheath)  that 
transitions  into  a  region  with  little  to  no  electric  field  (bulk  plasma).  These  conditions 
will  highlight  any  errors  in  code  implementation  very  quickly.  The  second  is  that  it  is 
well-documented  by  other  computational  references  (16:7-9;  36:4-5)  so  that  results  can 
easily  be  compared.  This  test  excludes  the  electron  energy  equation  assuming  constant 
transport  and  rate  coefficients. 

For  this  test,  an  Argon  plasma  is  modeled  at  a  pressure  of  100  torr.  The  transport 
and  rate  coefficients  are 


jue  =  0.3m2  tV ! s  =  10~3m2  IV I  s  V{  =  0.0 
D  =  0.3 m1  Is  D  =  10  4m2  Is  Pe  ,  =  0.0 

e  // 


(102) 


This  models  a  non-ionizing  plasma  with  a  characteristic  electron  temperature  held 
constant  at  \eV  and  an  ion  temperature  of  0. \eV  .  The  recombination  losses  are 
considered  to  be  negligible  for  this  problem. 
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Initially,  the  left  electrode  is  at  OV  and  the  right  side  of  the  domain  at  a  distance  of 
200  Debye  lengths  is  free  plasma.  The  Debye  length  XD0  is  based  on  conditions  in  the 
unperturbed  bulk  plasma  and  is  detennined  by 


where  Te  has  units  of  eV  .  The  electron  and  ion  densities  are  equivalent  and  spatially 

uniformat  ne  =n  =n0  =  1 .0  x  I  O'm  where  n0  is  the  initial  charge  density.  Because 

these  densities  are  equal  and  the  left  electrode  is  held  at  0V  ,  the  electric  potential  and 
field  are  initially  zero  as  well. 

At  time  t  =  Os ,  the  potential  of  the  left  electrode  is  reduced  to  -  50V  .  Formation 
of  the  sheath  is  initialized  as  electrons  are  repelled  and  ions  are  attracted  to  the  negatively 
biased  electrode.  Because  the  electrons  are  much  more  responsive,  they  move  quickly  to 
the  right  creating  a  large  charge  discrepancy  at  the  left  boundary.  As  time  progresses,  the 
charges  redistribute  in  such  a  way  as  to  try  to  neutralize  the  field  (16:7).  Figs.  21  and  22 
show  the  progression  of  the  charge  densities  and  the  electric  potential  respectively.  Note 
that  the  densities  are  shown  in  normalized  form  ns  /  n0  and  the  time  progression  is  given 

in  terms  of  t  /  tdi  where  tdi  is  the  dielectric  relaxation  time  in  terms  of  ions  only 


(104) 
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x!  A, do 


Figure  2 1 .  Transient  Sheath  -  Sheath  Progression 


Figure  22.  Transient  Sheath  -  Electric  Potential  Progression 


Within  several  tens  of  dual-particle  dielectric  relaxation  times  calculated  by  Eq. 
(70),  electrons  are  lost  to  the  bulk  plasma  at  the  right  edge  of  the  computational  domain 
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and  ions  are  pulled  from  the  bulk  plasma  and  lost  to  the  electrode  at  the  left.  The  highly 
mobile  electrons  become  distributed  in  such  a  way  that  cancels  out  the  field  for  the 
remainder  of  the  bulk  plasma  and  an  electric  potential  such  as  that  shown  for 
t  /  tdi  p  =  1 00  becomes  the  standard.  The  entire  process  continues  as  time  advances,  the 

sheath  moves  towards  the  bulk  plasma  and,  as  can  be  noted  from  Fig.  21,  the  sheath 
progression  slows. 

Qualitative  comparison  between  the  results  shown  in  Refs.  16,  36  and  Figs.  21 
and  22  shows  excellent  agreement.  This  similarity  demonstrates  the  accuracy  of  the 
semi-implicit  code  in  modeling  the  sheath  region  of  the  discharge.  While  the  transient 
sheath  test  case  excludes  the  energy  equation,  it  validates  Poisson’s  equation  coupled 
with  the  electron  and  ion  continuity  equations. 


Figure  23.  Ambipolar  Diffusion  (4:28) 


Ambipolar  Decay 

With  the  transient  sheath  model  validated,  the  next  test  case  to  be  analyzed  is 
ambipolar  decay.  Because  the  electrons  move  so  much  more  quickly  than  the  heavy  ions, 
an  initially  quasi-neutral  plasma  will  naturally  experience  a  charge  separation  (9:28). 
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When  there  is  sufficient  charge  density  to  generate  a  considerable  space  charge  as  they 
separate,  a  polarization  field  is  created  that  opposes  further  separation  (Fig.  23).  The 
charges  and  field  readjust  so  that  the  electrons  are  held  back  as  the  ions  are  pushed 
forward  -  they  can  only  diffuse  as  “a  team”. 

It  is  possible  to  derive  an  analytic  solution  to  both  the  particle  densities  and  the 
electric  field.  While  some  of  the  basics  will  be  covered  here,  Ref.  37  shows  greater  detail 
of  the  derivation  (131-133).  The  derivation  of  the  ambipolar  electric  field  begins  with  the 
assumption  of  quasi-neutrality:  ne  ~  n  ~  n  and  f  e  «  f  «  f  .  These  relationships  yield 

jupnE  -  DtVn  =  juenE  -  DeVn  .  (105) 


Solving  for  the  electric  field  gives  an  analytic  solution  of  the  form 

_  (Dp-D)s/n 

(Mp+Me)  n 


(106) 


The  particle  density  derivation  begins  by  substituting  this  electric  field  back  into 
the  ion  flux  equation  yields 


K--p.) 

(^+«) 


Vn  -  DVn  = 


/  MpDe  +  MeD,  ' 

Mp+Ee 


Vn . 


(107) 


Solving  the  electron  flux  equation  will  give  the  same  equation.  The  gradient  multiplier 
becomes  the  ambipolar  diffusion  coefficient: 
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MpPe  +  MeDp 
Mp  +  Me 


(108) 


Now  the  continuity  equation  can  be  written  as 


— +  V-(-D  Vn)=0.  (109) 

dt  V  ’ 


Because  each  of  the  coefficients  on  the  right  side  of  Eq.  (108)  is  spatially  uniform 
for  this  test  case,  Da  can  be  moved  through  the  first  gradient  operator  in  Eq.  (109). 
Assuming  that  the  solution  to  Eq.  (109)  can  be  separated  as 

n(x,t)  =  X{x)T{t)  (HO) 


then  this  equation  can  then  be  cast  into  the  form 


X 


dt 


=  DaT 


d2X 

dx1 


(HI) 


Dividing  Eq.  ( 1 1 1)  by  XT ,  the  variable  separation  is  complete.  Setting  each  equation 
equal  to  the  constant  -  Hr  and  solving  for  the  time  dependence  yields 


it)= 


n0  exp 


VU 


(112) 


The  solution  to  the  spatial  equation  is 


f 

=  sin 

V 


n 


\ 

) 


(113) 
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Since  the  density  is  modeled  at  zero  at  x  =  0  and  x  =  L  ,  this  means  that 


or  substituting  A  =  L  /  n  and  rearranging: 


r 


(114) 


(115) 


where  r  becomes  the  characteristic  time.  Putting  all  the  pieces  together,  the  combined 
solution  becomes 


i(x,t)  = 


n0  exp 


f  Da  ^ 

'  7Vp 

- ~t 

sin 

— 

{  A2  J 

UJ 

(116) 


For  this  test  case,  a  helium  plasma  is  modeled  at  1  torr  between  two  electrodes 
where  each  electrode  is  grounded.  The  initial  density  is  set  to  n0  =  1  x  10l6m  .  Since  the 
solution  will  be  in  the  form  of  a  sinusoid,  the  charge  densities  are  equal  and  initially  set  at 


ne{x)  =  n  (x)  =  n0Sin 


V  L  j 


(117) 


The  transport  parameters  are  held  constant  at 

/4  =  89.8m2  IV  Is  jup  =  0.8m2  IV Is  v,  =  0.0 
De  =  101.25m2  Is  Dp  =  (0.02 eV)jup  Pe  i  =  0.0 


(118) 


65 


modeling  a  non-ionizing  helium  plasma  with  a  characteristic  electron  temperature  of 
l.leV  and  an  ion  temperature  at  1/ 40 eh  .  Once  again,  the  energy  continuity  equation  is 
excluded  and  the  recombination  losses  considered  to  be  negligible  for  this  problem. 


Figure  24.  Ambipolar  Particle  Densities  -  Analytic  and  Computational 


x  (m) 


Figure  25.  Ambipolar  Electric  Field  -  Analytic  and  Computational 


66 


time(s) 

Figure  26.  Natural  Log  of  Peak  Density  -  Ambipolar  Decay 

At  time  t=0,  this  quasi-neutral  system  is  allowed  to  diffuse.  Results  were 
compared  after  one  characteristic  time  interval  r  .  At  first  glance,  the  results  did  not 
show  good  agreement  with  the  analytic  solutions.  As  Fig.  24  shows,  the  particle  density 
decays  more  quickly  than  the  analytic  solution  predicts.  Fig.  25  shows  that  the 
computational  electric  field  does  not  match  up  near  the  electrode  boundaries.  Lastly, 
since  the  decay  of  the  density  profile  is  governed  by  the  exponential  in  Eq.  (116),  the 
natural  log  of  the  density  at  one  particular  x-location  as  a  function  of  time  should  be 
linear  with  a  slope  of  -  Da  /  A2 .  Fig.  26  shows  the  natural  log  of  the  peak  density  as  a 

function  of  time  along  with  the  least  squares  linear  fit.  While  the  computational  results 
are  close  to  linear,  there  is  a  slight  curve  to  the  results.  The  slope  of  the  least  squares  line 
is  significantly  more  negative  than  the  expected  value  of  -5636.1 . 
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Upon  further  analysis,  it  was  determined  that  the  discrepancies  came  from  the 
existence  of  sheath  regions  near  both  electrodes.  This  diagnosis  came  from  the  fact  that 
as  the  initial  density  was  increased,  the  slope  of  the  natural  log  plot  became  closer  to  the 
analytic  value.  Increasing  the  initial  density  decreases  the  size  of  the  sheath  region  and 
therefore  decreased  the  size  of  the  problem  area. 

The  sheath  is  a  problem  area  because  as  electrons  and  ions  near  the  edge  of  the 
discharge,  the  charge  density  greatly  diminishes  and  free  diffusion,  as  opposed  to 
ambipolar,  now  applies  (9:28).  The  field  created  by  the  space  charge  separation  is  no 
longer  large  enough  to  keep  the  charged  particles  together  and  the  electrons  quickly  leave 
the  ions  far  behind.  Quasi-neutrality  no  longer  exists.  Since  the  ambipolar  analytic 
results  depend  on  the  maintenance  of  quasi-neutrality,  the  equations  derived  above  do  not 
apply  in  this  region. 

Further  proof  of  the  sheath’s  accelerating  effect  on  the  decay  rate  came  with  the 
quantification  of  the  sheath  size.  The  density  decay  rate,  assumed  to  be  a  constant,  is 
given  as 


1  Dg  Da^ 

r  A2  L2 


(119) 


This  equation  holds  if  the  particle  densities  between  x  =  0  and  x  =  L  remain  quasi¬ 
neutral.  As  time  progresses  the  region  of  quasi-neutrality  decreases.  Using  the  Debye 
length  as  the  characteristic  length  gives  the  time-dependent  equation  for  the  effective  A : 

A^=L~2a-AD{t)  (120 

n 
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where  a  represents  the  effective  sheath  thickness  and  the  time-dependent  Debye  length 
relates  back  to  the  initial  Debye  length  from  Eq.  (103)  as 

-U'Mm  J^j-  (120 

The  time-dependent  peak  density  was  used  as  the  n(t)  in  Eq.  (121).  The 
computational  decay  rate  was  then  fit  using  Eq.  (120)  in  order  to  solve  for  a  .  A  sheath 
thickness  equal  to  7.66  gave  results  that  matched  the  peak  density  decrease  with  less  than 
1%  error.  Figs.  27  and  28  show  these  results. 


time  (s) 


Figure  27.  Time-Dependent  Peak  Densities 
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time  (s) 


Figure  28.  Percent  Error  for  Peak  Density  Decay 

This  analysis  shows  excellent  agreement  between  the  computational  and  analytic 
results.  The  ambipolar  test  case  was  performed  using  both  the  semi-implicit  and  fully 
implicit  numerical  codes.  While  the  semi-implicit  results  are  shown  above,  the  fully 
implicit  showed  a  peak  density  difference  of  1.85%  after  the  characteristic  time  and 
showed  excellent  agreement  in  all  regions.  This  test  validated  both  codes  as  to  the 
implementation  of  Poisson’s  equation  coupled  to  the  electron  and  ion  continuity 
equations.  Now  that  the  first  three  equations  in  the  system  have  been  solved  and 
validated,  the  energy  equation  will  be  tested. 

Radio-Frequency  Glow  Discharge 

For  the  final  test  of  the  semi-implicit  code,  a  radio-frequency  source  was 
implemented  along  with  the  energy  equation  and  source/loss  terms.  Similar  to  the 
transient  sheath,  this  problem  is  an  excellent  initial  test  of  the  energy  equation  and 
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source/loss  implementation  for  two  reasons.  First,  the  alternating  sheath  will  be  an 
excellent  test  of  the  code’s  ability  to  account  for  change  in  regions  of  strong  fields. 
Second,  this  test  is  also  well  documented  by  other  resources  (16:10-12;  18;  38:2785- 
2786)  and  can  therefore  be  qualitatively  compared. 

For  this  test  case,  a  helium  plasma  is  modeled  at  1  torr  between  two  electrodes 
separated  a  distance  of  0.04m.  The  potential  at  the  left  electrode  oscillates  as  a  sinusoid 
with  a  frequency  of  10  MHz  and  an  amplitude  of  500V 

<pL  (t)  =  500sin(2^/i)  (122) 

and  the  right  electrode  is  grounded.  The  transport  and  source/loss  coefficients  are 

Me  =  Me  k]  Mp  =  Mp[e  IP]  v,  =  v[ue] 

De=De[ue]  Dp={0.02eV)<up  =  k-Om'3 ' 

The  electron  mobility  and  diffusion  and  the  ion  mobility  coefficients  are  obtained  by 
piecewise  fits  to  the  curves  shown  in  Figs.  8,  9,  10  and  1 1  (Appendix  C). 

For  this  test,  the  initial  densities  were  again  equal  and  assumed  to  have  a 
sinusoidal  spatial  distribution  as  shown  in  Eq.  (117)  with  n0  =  5  x  1015  m  .  At  t=0,  the 

AC  voltage  source  was  turned  on  and  the  discharge  was  allowed  to  run  until  significant 
changes  in  the  peak  density  terminate.  This  takes  approximately  400  cycles  for  this 
simulation.  The  results  were  recorded  after  the  discharge  had  been  in  steady-state  for 
more  than  100  cycles.  A  qualitative  comparison  of  the  particle  densities,  electric  field 
and  current  densities  are  then  made  to  the  results  shown  in  Refs.  16,  18  and  38. 
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Figure  29. 


Semi-Implicit  RF  Densities  and  Electric  Fields 


Figure  30.  Results  from  Hilbun  -  RF  Densities  and  Electric  Fields  (16:12) 


72 


Densityi^O11' m  3).3,  w  Density(l0'6 m 


Because  the  density  results  from  each  of  the  references  are  qualitatively  similar, 
only  the  results  from  Hilbun  in  Ref.  16  are  shown.  Fig.  29  shows  the  computational 
results  obtained  through  this  numerical  simulation  and  Fig.  30  shows  the  results 
presented  by  Hilbun.  The  result  profiles  show  excellent  agreement.  Note  that  the  results 
presented  are  a  quarter  cycle  out  of  phase  due  to  a  phase  shift  in  the  driving  frequency 
and  that  this  code  reports  the  values  in  terms  of  SI  units  whereas  the  Hilbun  reference 
gives  them  in  CGS. 

For  this  test  once  in  steady-state,  the  ion  densities  change  very  little  over  one 
voltage  cycle.  The  mobile  electrons,  however,  move  quickly  in  and  out  of  the  sheath 
regions  of  the  discharge.  As  it  cycles,  they  are  both  repelled  by  the  cycle-dependent 
anode  and  attracted  to  the  cycle-dependent  cathode.  At  the  quarter  and  three-quarter 
cycle  points,  the  electrons  move  out  of  the  cycling  cathode  sheath  and  the  positive  space 
charge  formed  creates  a  significant  electric  field.  At  the  zero  and  half  cycle  points 
(measured  by  this  project)  when  the  potential  is  zero  on  both  electrodes,  the  electrons  are 
evenly  distributed  with  a  small  positive  space  charge  creating  a  diminished  electric  field 
on  each  side  of  the  discharge  regime. 

For  the  current  density  comparison,  the  electron,  ion  and  displacement  current 
densities  were  tracked  over  one  cycle  at  the  left  (driven)  electrode.  These  current 
densities  are  found  with  the  equations 


J 


e,U2 


=  qe  • 


,  ne,lVe,th,H2 

V  4 


(124) 
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Jp,i/2  ~  Qpnp,\MP, \E\n 


(125) 


and 


J d,m  ~  £o ' 


dE 


1/2 


dt 


(126) 


where  q  is  the  signed  electron  charge.  As  can  be  seen  from  Fig.  31,  the  displacement 
current  density  at  this  location  accounts  for  most  of  the  total  current  density  for  all  phases 
of  the  cycle. 


Figure  3 1 .  RF  Current  Densities  at  Left  Electrode 

Comparison  of  current  densities  between  this  project  and  the  other  references 
once  again  showed  excellent  agreement.  There  were  some  differences  between  the 
shapes  of  each  individual  contribution  (Figs.  3 1  and  32).  The  fact  that  the  electron  flux 
for  this  code  is  modeled  as  thermal  at  the  boundaries  accounts  for  one  of  the  primary 
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differences.  The  fact  that  the  total  current  density  must  be  constant  and  should  be  cyclic, 
phase-shifted  from  the  voltage  (38:2786)  means  that  a  change  to  the  electron  current 
density  must  change  the  displacement  and  ion  current  densities  as  well.  The  total  current 
density,  however,  was  nearly  identical  in  both  shape  and  magnitude  for  each  of  the 
comparative  references. 


F  =  10MHz 


Figure  32.  Results  from  Hilbun  -  RF  Current  Densities  (16:13) 


Besides  small  differences  that  could  be  accounted  for  with  different  or  absent 
treatment  of  the  energy  equation  as  well  as  the  thermal  boundary  flux  conditions  for 
electrons,  this  code  showed  excellent  agreement  to  the  references.  With  this  test 
completed,  the  full  semi-implicit  model  has  been  validated.  The  DBD  effect  on 
discharge  current  can  now  be  accurately  characterized. 
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The  DBD  Streamer 


The  discharge  behavior  changes  when  a  dielectric  barrier  covers  one  of  the 
electrodes.  As  a  way  to  study  this  configuration,  this  project  considers  streamer  (29:93) 
or  streamer-regime  filament  (39:7569)  fonnation  in  a  DBD  configuration.  A  streamer  is 
essentially  an  ionization  wave  crossing  the  discharge  gap.  These  waves  will  typically 
form  when  a  voltage  that  is  significantly  greater  than  the  breakdown  voltage  is  applied  to 
a  discharge  (39:7569).  While  much  of  the  computational  and  experimental  DBD  studies 
are  conducted  using  an  AC  source,  this  project  will  only  address  the  DC  simulation  of 
this  process.  Because  the  current  spike  occurs  so  quickly  in  an  alternating  current  DBD 
that  the  driving  voltage  changes  very  little  during  the  process,  a  DC  source  is  a 
reasonable  approximation. 

For  this  simulation,  the  semi-implicit  code  was  used  to  model  a  helium  plasma  at 
1 0  torr.  Initially,  the  charged  particle  densities  are  set  to  a  very  low  value 
( n0  =  5  x  109  m  )  with  an  ionization  fraction  ~10"13.  The  initial  charge  density 

configuration  resembles  that  of  a  discharge  with  a  cathode  sheath  already  formed  (Fig. 
34a).  Because  the  assumed  low  value  of  the  initial  charge  densities,  the  electric  field  is 
essentially  constant  throughout  the  computational  domain  at  the  initiation  of  the 
simulation.  The  transport  and  rate  coefficients  used  are  the  same  as  those  for  the  RF  test 
case  given  in  Eq.  (123). 

An  additional  loss  term  is  added  to  the  system  of  equations  in  order  to  better 
model  an  actual  discharge.  This  term  considers  a  diffusion  loss  to  the  walls  of  the 
discharge  tube.  This  radial  diffusion  is  considered  to  be  ambipolar  and  the  loss  is  taken 
as 
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(127) 


Z  —£*-n 

^ Da  ^2  Hs 


where  the  ‘s’  makes  the  loss  specific  to  the  density  of  each  particle  type.  Even  though 
this  loss  is  modeled  as  ambipolar,  the  charged  particle  densities  will  be  different  in 
certain  regions  of  the  discharge  and  this  general  equation  takes  that  into  account.  Since 
this  model  only  considers  radial  losses,  A  =  A/2.4  (9:67)  where  R  represents  the  radius 
of  the  discharge  tube  at  0.02m. 


I  Vgap  I 


Figure  33.  Dielectric  Circuit  Configuration  for  Streamer  Simulation 


For  the  voltage  configuration  of  this  test  case,  a  circuit  is  constructed  to  resemble 
that  shown  in  Fig..  The  applied  system  voltage  is  -2.5kV  and  the  ballast  resistance  is  set 
to  12.5kf2.  The  effective  capacitance  of  the  dielectric  barrier  is  assumed  to  have  a 
negligible  effect  on  the  circuit  behavior  and  is  therefore  ignored.  The  effect  of  the 
dielectric  barrier  is  modeled  as  described  by  Eq.  (58). 

The  potential  applied  to  the  left  electrode  becomes 

Vgap  =V applied  _  ^ ballast  (128) 
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where  /  represents  the  current  determined  by 

/  =  (129) 

A  is  the  area  of  the  left  electrode  assumed  to  be  (0.05)';r .  The  current  densities  are 
determined  by  Eqs.  (124)  and  (125)  and  is  calculated  at  x  =  0.02m  .  The  displacement 
current  is  excluded  because  of  difficulties  incorporating  this  contribution  into  the 
feedback  system  implemented  for  this  simulation  and  is  recognized  to  be  a  source  of 
some  error  in  the  results  presented.  The  general  features  of  the  discharge  evolution, 
however,  are  consistent  with  experimental  observations  and  earlier  simulations  in  a  bare 
electrode  configuration  (16:15-18).  The  right  electrode  is  grounded. 


(c)  t  =  48ns 


(d) t  =  1 14ns 
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Electric  Field  (kV/m)  Electric  Field  (kV/m) 


Figure  34.  DBD  Streamer  Evolution  -  (a)  Initial  distribution  (b)  Homogenous  avalanche 
(c)  Space-charge  dominated  avalanche  (d)  Cathode-directed  streamer 

The  results  of  this  simulation  show  the  evolution  of  the  streamer  and  the  effect  of 
the  dielectric  barrier  on  this  evolution  (Fig.  34a-d)  .  Because  the  simulation  begins  with  a 
plateau  of  charge  across  the  right  side  of  the  domain,  the  density  growth  and  sheath 
progression  to  the  right  represent  the  electron  avalanche  associated  with  breakdown  in  the 
gap  (Fig.  34b).  The  wave  of  ionization  is  directed  toward  the  anode  (buried  electrode)  in 
this  phase.  As  the  electrons  are  driven  to  the  right  by  the  negative  (left-directed)  electric 
field,  the  ionization  cascade  creates  a  trail  of  ions  that  can  not  quite  keep  up  with  the 
electron  progression.  This  phase  of  the  discharge  cycle  is  called  the  homogenous 
avalanche  (16:17)  because  the  density  growth  is  governed  by  the  electron  advance  in  a 
nearly  uniform  electric  field. 

The  next  phase  in  the  streamer  evolution  occurs  as  the  charge  densities  and  space 
charge  separation  become  significant  enough  to  alter  the  electric  field  (Fig.  34c).  The 
charges  redistribute  themselves  so  as  to  reduce  the  electric  field  in  the  avalanche 
(maximum  charge)  region  of  the  discharge.  As  these  charges  diminish  the  electric  field 
within  the  avalanche,  the  field  to  the  left  and  at  the  dielectric  becomes  even  stronger. 
During  this  phase,  a  shift  in  the  ionization  occurs.  Originally,  the  greatest  ionization 
occurred  within  the  avalanche  as  it  progressed  from  the  cathode  to  the  anode.  Now,  more 
ionization  occurs  to  the  left  of  the  avalanche  where  the  field  has  increased  than  in  the 
avalanche  itself  where  the  field  has  weakened.  This  is  the  space-charge  dominated 
avalanche  phase  (16: 17). 


79 


The  final  phase  is  associated  with  a  change  in  the  direction  of  the  wave  now 
streaming  towards  the  cathode  (Fig.  34d).  As  the  field  strength  to  the  left  of  the 
avalanche  continues  to  grow,  the  ionization  in  this  region  remains  significant.  Because 
the  field  is  still  negative,  it  pulls  the  new  ions  to  the  left  as  it  pushes  the  new  electrons  to 
the  right.  Accordingly,  this  new  ionization  wave,  labeled  the  cathode-directed  streamer 
(16:17),  has  a  net  positive  charge  at  the  streamer  head.  This  positive  space  charge  acts  as 
an  effective  anode  and  as  the  streamer  continues  to  approach  the  cathode,  the  electric 
field  becomes  increasingly  negative  in  the  interim  region. 

During  the  first  two  phases  of  the  discharge  (Fig.  34b  and  34c),  the  avalanche 
progression  is  very  similar  to  that  of  a  non-barrier  microdischarge  (16:16).  The  surface 
charge  on  the  dielectric  increases  during  the  avalanche,  however,  and  alters  the  discharge 
behavior  for  the  final  phase  of  the  DBD.  During  this  phase,  the  charge  accumulation 
plays  an  important  role. 

As  the  ionization  wave  from  an  anode-directed  avalanche  transitions  to  a  cathode- 
directed  streamer,  the  negative  surface  charge  that  has  accumulated  on  the  dielectric 
reaches  a  saturation  value.  This  surface  charge  is  shown  in  Fig.  (34)  as  the  right-most 
point  in  each  of  the  density  plots  and  should  be  distinguished  from  the  volume  charge 
density  shown  for  the  remainder  of  the  domain.  The  particles  between  the  dielectric 
surface  and  the  streamer  head  begin  to  redistribute  such  that  the  space  charge  diminishes 
and  the  field  in  this  region  continues  to  diminish  as  well. 

This  field  reduction  becomes  important  when  considering  the  discharge  circuit 
described  by  Eq.  (128).  As  the  avalanche  progresses  to  the  right,  the  current  flowing  to 
the  left  electrode  diminishes  the  gap  potential.  As  the  charge  density  grows  increasingly 
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negative  on  the  dielectric  surface,  the  gap  potential  is  reduced  even  further.  A  sharp  rise 
in  the  current  associated  with  the  cathode-directed  streamer  again  reduces  the  gap 
potential.  This  reduction  in  the  gap  potential  will  progress  due  to  the  accumulation  of 
surface  charge  until  the  voltage  falls  below  the  self-maintenance  value  and  the  discharge 
will  eventually  extinguish.  In  this  relatively  low  pressure  discharge  where  recombination 
is  not  as  large,  this  progression  may  include  a  temporary  transition  to  a  glow  discharge 
before  this  charge  build-up  on  the  dielectric  becomes  significant  enough  to  extinguish  the 
discharge. 

The  importance  of  streamer  or  filament  formation  in  the  DBD  becomes  evident 
when  comparing  the  right-ward  and  left-ward  movement  of  the  charged  particles. 

Because  the  streamer  carries  both  ions  and  electrons  in  its  ionization  wave,  the 
asymmetry  associated  with  the  streamer  is  a  very  interesting  characteristic.  There  are  a 
significantly  larger  number  of  ions  involved  with  the  cathode-directed  streamer  than  there 
are  in  the  anode-directed  avalanche.  If  it  is  assumed  that  ion-neutral  collisions  dominate 
the  momentum  transfer  to  the  flow  (16:18)  then  a  source  of  the  observed  asymmetry  in 
flow  control  could  be  explained,  at  least  in  part,  by  streamer  propagation. 

The  true  determination  of  the  source  of  the  added  flow  momentum  has  yet  to  be 
definitively  determined.  There  is  a  great  deal  more  research  that  remains  to  be  done  on 
this  subject.  Additionally,  there  are  several  model  assumptions  made  in  this  project  that 
deserve  more  investigation.  Each  of  these  subjects  will  be  covered  in  the  final  chapter. 
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V.  DBD  Modeling  -  Accomplishments  and  the  Next  Round 


This  numerical  study  examined  plasma  dynamics  in  a  dielectric-barrier  discharge 
configuration.  While  this  investigation  took  one  small  step  towards  a  better 
understanding  of  DBD  operation,  there  still  remains  a  great  deal  more  to  be  studied  and 
characterized.  Some  of  the  lessons  learned  as  well  as  ideas  for  continuing  investigation 
will  be  discussed. 

Looking  Back 

This  research  achieved  far  more  than  a  simple  DBD  streamer  characterization. 

The  intennediate  results  of  this  research  were  validated  against  previously  published 
computational  models.  An  accurate  one-dimensional  code  was  implemented,  tested  and 
validated.  An  alternative  numerical  approach  was  also  implemented  and  partially 
validated.  This  fully  implicit  method  improved  upon  the  semi-implicit  formalism  with  an 
increased  time  step  and  the  promise  of  increased  accuracy.  The  validation  procedure  and 
observed  consistency  of  the  two  numerical  approaches  ensured  that  the  DBD  simulation 
yielded  accurate  numerical  results.  This  simulation  extended  current  numerical  work  in 
characterizing  the  streamer  cycle  when  a  dielectric  barrier  blocks  one  electrode. 

Accurate  1-D  Model 

Before  the  computational  models  could  be  used  to  investigate  DBD 
characteristics,  they  first  had  to  be  validated.  While  there  is  much  documentation  as  to 
the  numerical  intricacies  involved  in  developing  a  plasma  discharge  model,  this  specific 
implementation  had  to  be  rigorously  tested  and  the  results  compared  to  either  previously 
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validated  results  or  analytic  solutions.  The  full  implementation  of  the  semi-implicit  code 
was  tested  and  validated  through  the  analytic  Poisson,  transient  sheath,  ambipolar  decay 
and  RF  simulations. 

The  Poisson  check  showed  excellent  agreement  with  the  analytic  results  and 
demonstrated  that  the  solution  to  Poisson’s  equation  was  accurate  for  both  the  semi- 
implicit  and  fully  implicit  models.  The  transient  sheath  test  case  agreed  with  two 
previously  validated  models  and  showed  that  the  coupling  between  the  electron/ion 
continuity  equations  and  Poisson’s  had  been  correctly  implemented.  It  also  established 
that  the  semi-implicit  code  yielded  an  accurate  model  of  the  discharge  sheath.  The 
ambipolar  test  case  showed  excellent  agreement  with  the  analytic  solutions  for  the 
electric  field  and  charge  density  profile  once  the  sheath  region  was  taken  into  account. 
The  sheath  analysis  showed  agreement  to  within  1%  relative  error  for  every  time  step  of 
the  decay  and  validated  the  implementation  of  the  ion/electron  continuity  equations 
coupled  with  Poisson’s  for  both  the  semi-implicit  and  fully  implicit  models. 

Fully  Implicit  Advantages 

Although  there  was  not  enough  time  to  completely  implement  and  test  the  fully 
implicit  code,  the  benefits  of  this  solution  technique  became  obvious  in  the  ambipolar 
diffusion  test.  Because  there  is  no  lag  between  Poisson’s  solution  and  the  particle/energy 
continuity  calculations,  a  larger  time  step  can  be  taken.  For  the  ambipolar  test  case,  this 
time  step  was  three  to  four  orders  of  magnitude  greater  than  that  of  the  semi-implicit 
method.  Additionally,  the  total  execution  time  was  reduced  by  a  factor  of  four. 


83 


DBD  Streamer  Characterization 


The  final  phase  to  this  project  simulated  the  streamer  formation  during  a  DBD 
discharge  initiation.  This  characterization  detailed  three  different  phases  of  the  DBD 
cycle.  The  first  two  phases  of  the  DBD  were  the  anode-directed  homogenous  avalanche 
and  space-charge  dominated  avalanche.  The  third  phase  was  the  cathode-directed 
streamer.  This  simulation  confirmed  the  asymmetry  between  the  anode-directed 
avalanche  and  the  cathode-directed  streamer.  A  comparison  to  experimental  and  non¬ 
barrier  computational  results  showed  general  agreement. 

Looking  Forward 

Fully  Implicit  Development 

Due  to  time  constraints,  a  full  implementation  for  the  simultaneous  solution  could 
not  be  accomplished.  Implementation  of  the  energy  equation  was  problematic.  In 
regions  where  the  electric  field  was  strong  and  the  electron  densities  were  small,  negative 
energies  were  encountered  in  the  Newton-Raphson  iteration  process.  Additionally,  once 
the  energy  equation  was  implemented,  the  time  steps  required  to  maintain  stability 
dropped  below  those  used  for  the  semi-implicit  solution. 

Boeuf  introduces  a  time-dependent  Poisson’s  equation  in  Ref.  25.  In  that 
particular  study,  the  equation  is  used  in  a  semi-implicit  numerical  formalism.  However, 
Boeuf  reports  that  this  implementation  added  enough  stability  to  the  method  that  time 
steps  as  large  as  fifty  times  the  dielectric  relaxation  time  could  be  taken.  This  added 
stability  would  certainly  transfer  to  the  fully  implicit  model  as  well. 
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Mean  Velocity  Flux 

Hammond  proposes  a  novel  discretization  technique  in  his  paper  addressing  RF 
plasma  discharge  simulations  (Ref.  40).  This  mean  velocity  flux  discretization  is  a 
method  that  could  be  employed  to  replace  the  Scharfetter-Gummel  discretization.  This 
method  promises  increased  accuracy  of  the  solution  and  a  decrease  in  the  computational 
cost  of  the  flux  calculations.  This  method  of  flux  discretization  certainly  deserves  more 
exploration. 

Continued  Dielectric  Studies 

There  still  remain  a  large  number  of  topics  to  be  considered  pertaining  to  the 
dielectric  barrier  discharge.  Studies  involving  AC  voltage  sources,  differing 
dielectric/charge-interaction  models  as  well  as  moving  to  two  or  three  dimensions  and 
various  electrode  configurations  name  just  a  few  of  the  possibilities.  Something  that 
becomes  increasingly  important  as  the  investigation  moves  to  multiple  dimensions  is  the 
Poisson  solver.  If  the  asymmetric  electrode  configuration  shown  in  Fig.  7  is  modeled,  the 
boundary  conditions  for  Poisson’s  equation  get  far  more  complicated. 

In  order  to  pin-point  the  source  of  the  observed  flow  control,  many  of  the  subjects 
listed  above  will  need  to  be  addressed.  The  future  goal  for  the  continuation  of  this 
project  is  to  create  a  computational  tool  that  can  be  used  to  simulate  a  plasma  discharge 
in  all  configurations.  This  tool  would  include  the  ability  to  add  the  dielectric  barrier  to 
any  number  of  electrodes,  model  in  one  or  two  dimensions  with  a  variety  of  electrode 
geometries.  This  tool  could  be  used  to  study  the  effects  of  the  barrier  on  the  discharges. 

In  the  process,  it  could  possibly  pinpoint  the  mechanism  that  drives  DBD  flow  control. 
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Depending  upon  what  this  mechanism  is,  the  possibilities  for  actuator  technologies  could 
be  endless. 
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Appendix  A.  Model  Assumptions 


Many  of  the  model  assumptions  match  those  of  Boeuf  and  Pitchford  (14: 1377).  The 
general  model  assumptions  for  this  project  are: 

-  Neglect  interactions  between  the  charged  particles  and  neutral  atoms  in  excited 
states 

-  Only  electron  ionization  from  the  ground  state  is  considered 

-  Electron-electron  collisions  are  not  taken  into  account 

-  Ion  inertia  is  neglected 

-  Ion  distribution  is  assumed  Maxwellian 

-  Electron  inertia  and  energy  gradient  terms  are  neglected  in  the  electron 
momentum  equation 

-  Pressure  tensor  is  assumed  isotropic  and  diagonal 

-  Electron  drift  energy  is  considered  negligible  with  respect  to  the  electron 
thennal  energy 

-  Heat  flux  is  proportional  to  the  electron  temperature  gradient 

-  Mean  electron-neutral  collision  rates  are  proportional  only  to  the  electron  mean 
energy 

-  The  electron  diffusion  term  is  essentially  constant  over  small  spatial  regions  and 
can  be  pulled  through  the  gradient  operator  as  a  constant  for  the  drift  portion  of 
the  energy  density  flux  equation 
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Appendix  B.  Newton-Raphson  Matrix 


Original  PHI  equation: 


PHI[cp,  ne,  np\ 


f  ( 1jW+1  \  r\(  k+\,m+\  \  .  (  k+l,m+l\\ 

I  Wm  )~2Wi  )+\<Pi-i  j 


(Ar)2 


+  «-(n 


k+\,m+l 

p,i 


-n 


k+\,m+\ 


)  =  0 


130) 


Newton-Raphson  equation: 


PHI\<p,ne,np]  = 


r  1  ' 

(Ax)2 


k,m+ 1 


c*  &,m+l  H  o  k,m+l  ,  H  c*  k,m+\ 

°<Pi - ™e\  +—5npi 


(A*) 

(^r+1)-2(^-m+1)+(^r+1)' 

(Ax)2 


+— «r1-<r1)=o 


(131) 


88 


Original  NP  equation: 
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Newton-Raphson  equation: 
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Original  EN  equation  (not  yet  validated): 
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Newton-Raphson  equation: 
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Coefficients  for  At  -  first  index  is  row,  second  index  is  column  (all  fourth  columns  are 
not  yet  validated): 
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Coefficients  for  B, : 
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Coefficients  for  C, : 
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Values  for  6,  vector  (f>4  not  yet  validated): 
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Appendix  C.  Transport/Rate  Coefficients  Fit 

BOLSIG  fit  for  the  electron  mobility  (17): 

(i ue  <  0.543)  — »  nep  =  1 .6e6 

(0.543  <ue  <  1.6)— >  /Jep  =  -6.0830  le5  -ue  +  1.81062e6 
(1.6  <ue  <  5.207) ->  pep  =  (-23. 87736  +922.45227)3 

(5.207  <  ue  <  7.087)  pep  =  (3.487968  -ue  +785.6738)2 
(7.087  <  ue  <  10.08)  -*  pep  =  (8.084379  •  ue  +  752.824841)2 
(10.08  <  ue  <  21.85)->  pep  =  16794.7345 -ue  +5.3014e5 
(2 1 .85  <  ue  <  53.68)  ->  pep  =  6252. 12852  •  ue  +  7.6784e5 
(53.68  <  ue )  — >  pep  =  \  .45550e5  •  ln(nj+ 5. 17329e5 


BOLSIG  fit  for  the  electron  diffusion  (17): 

(ue  <  6.75)  Dep  =  3.65887e5  •  ue  +  4.63643e5 

(6.75  <  ue  <  10.0 )^Dep  =  -4.86929e5  -ue  -3.5917e5 

(lO.O  <ue<  25.0)  — »  Dep  =  8380.8  •  {ue)2  +  4.50769e5  •  ue  -  8.552 12e5 

(25.0  <  ue  <  100.0)  ->  Dep  =  8.70618e5  •  ue  -  7.0e6 

(100.0  <ue<  200.0)  ->  Dl,P  =  9.48006e5  •  ue  - 1 .4<?7 

(200.0  <ue<  385.0)  ->  Dep  =  \.05e6  ■  ue  -  3.0e7 

(385.0  <  «e)  — >  Dep  =  3.7eS 


Gas  Parameter  File  (SIGLO-RF)  fit  for  the  ion  mobility  (18): 

(E/P  <  10.0) -+MpP  =  -73.967  -(E/P)  +  8725.8 
(E/P>  10.0)  ->  npP  =  exp(9.9  )-(e/pY035226 


BOLSIG  fit  for  ionization  rate  (17): 
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( ue  <  2.53)  —>vi!p  =  0.0 

(2.53  <ue<  2.88) -+vt/p  =  0.02 

(2.88  <  ue  <  3.26)  —>vi/p  =  0.33 

(3.26  <ue<  7.34)  — >  v.  / p  =  1.3e3 -exp(2.3112- we) 

(7.34  <  M<r  <  8.45)  -+v,/p  =  3.47826e  -11-  (ue  )17'326 
(8.45  <  ue  <  9.47)  -*vi/p  =  8.57143e  -  6  •  (ue  J1 1,525 
(9.47  <  we  <  12.4) ->  /  p  =  0.2 138  .(i/J7- 0431 

(12.4  <  <  49.92)  — >  vi  I  p  =  7.  Ie6  •  Me  -  7.7e7 

(49.92  <ue)^vi/p  =  2.4eS  ■  ln(«e)-6.6e8 
v,[.S7]  =  (^ //;)/; 

BOLSIG  fit  for  ionization  rate  (17): 

(ue  =  0)  ->  k,  =0.0 

(«e  <5.21)-^  =  2.0e- 12  •  (ue)2  +  1.3e-l  1  •  we  -  2.6e- 12 

(5.21  <ue  <7.09)^kL  =4.37e-ll-(«e)3  -7.25e-10-(nJ2  +4.05e-9-«e -7.46e-9 

(7.09  <  we  <  9.8)  ->  kL  =  8.0e- 10  ■  (ue)2  +  1.1255e-8  •  ue  -  6.6e-  8 

(9.08  <ue)^>kL  =  4.23<?  -  7 

(8.45  <  <  9.47)  ->  v,.  /  =  8.57143e  -  6  •  (ue )n'525 

(9.47  <  ue  <  12.4)  -^vi/p  =  0.2138  •  (u  e)7°m 

(12.4  <  ue  <  49.92)  ->  vt/p  =  7 .  Ie6  ■  ue  - 1.1  el 

(49.92  <  «e)— »  ia  / p  =  2Ae8-ln(ue)-6.6e& 
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